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ABSTRACT 


A computer model for axially symmetric dielectric radomes in the near field of 
a circular aperture has been developed. This model, based on a method of 
moments solution for bodies of revolution, is used to determine the electric field 
pattern of a linearly polarized circular aperture radiating through a dielectric 
radome. The program is written in FORTRAN and can accommodate rotationally 
symmetric radome shapes. Arbitrary illumination distributions for the aperture can 
be specified, and phased scanning of the aperture in both azimuth and elevation 
is allowed. Computational run time has been reduced by taking advantage of 
mode symmetry. A separate program has been developed to compute gain, and 
together these programs are used to determine the pattern degradation and gain 


loss due to the presence of a radome in the near field. 
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I. INTRODUCTION 


A radome is a structure designed to protect an antenna in its operating 
environment. Aircraft radomes must protect radar antennas from high 
aerodynamic forces while minimizing drag, weight, and interference on its 
electrical operation. The presence of a radome can affect the gain, beamwidth, 
sidelobe level, and direction of boresight, as well as change the VSWR and 
antenna noise temperature. This thesis involves the continued development of a 
computer code that predicts antenna radiation patterns looking through 
aerodynamic radomes. 

The ability to accurately predict the effect of a radome on the operation of an 
antenna early in the design process is desirable. Typically, designers of the 
earliest radomes were concerned more with structural considerations than 
electrical performance. However, after the near failure of the U.S. Army's 
“Pathfinder’ program of radar-equipped B-24 "Liberator" bombers due to poor 
electrical peiformance of their radomes, it became clear that," There was a serious 
need for developing sound theoretical procedure for the electrical design of 
radomes." [Ref. 1] For most practical radomes the surface geometry is too 
complex to yield closed-form mathematical expressions describing their 
performance. It is possible to predict radome performance using numerical 
solutions of Maxwell's equations. The program described in this thesis is based 
on a numerical technique called the method of moments (MM). The particular 
form of the solution used here is based on a formulation by Mautz and Harrington 
[Ref. 7] for bodies of revolution. Francis [Ref. 2] developed a computer code 


LDBORMM.F that specialized the Mautz and Harrington solution to ogive- 


shaped radomes. There are no restrictions on the location of the radome with 
respect to the antenna, and therefore near-field radomes can be analyzed. 

The goal of this thesis was to extend the capabilities of this program. 
Significant improvements include: taking advantage of symmetry to reduce 
computation time, applying realistic radar antenna patterns for comparison to 


experimental results, and the computation of gain loss due to the radome. 


A. DESCRIPTION OF THE PROBLEM 

The basic system analyzed is a radome consisting of a body of revolution 
(BOR) located coaxial with a circular antenna, as in Figure 1.1. The BOR 
constraint is satisfied by many aerodynamic radomes and greatly simplifies the 


MM solution , which will be described in more detail in following chapters. The 


Axis of 
Symmetry 


Figure 1.1 Radome and Antenna 


ogive shape of the radome in Figure 1.1 is the primary one of interest, but the 
program can also model other bodies of revolution including: a cone, disk, 
hemisphere, and parabola, as shown in Figure 1.2. These were primarily added for 


validation purposes. 


There are no constraints on the size, shape and location of these objects 


provided they are coaxial with the antenna. However, because MM requires the 


Figure 1.2 Optional Radome Shapes 


solution of a matrix equation, computer run time and memory will limit the size of 
the bodies that can be analyzed. A unique feature of this program is that it takes 
into consideration the fact that the radome may be in the near field of the antenna 
as in Figure 1.3. The antenna aperture used in the program is circular, but the 
program can be easily altered to accommodate other shapes. There are no 
constraints on the aperture size, scan angle, or polarization. 

The remainder of this thesis is divided into three major sections. Chapter II 
covers theoretical development of the problem and method of solution. The MM 
solution is described in Chapter II along with a discussion of mode symmetry and 
how it can be used to advantage in obtaining the radome current. Chapter II also 
describes how the total field is calculated from the currents. Chapter III covers 


the programs used for generating patterns based on the methods of 


Chapter II. The first part describes the program RADOME.F which calculates the 


current and electric field. The second part describes the program GAIN.F which 





Figure 1.3 Radome Located in Near Field of Antenna 


calculates the gain of the system. The final chapter discusses the results and 


presents several conclusions based on the calculated radome patterns. 


Il. ANALYSIS 


A. GENERAL DESCRIPTION OF ANALYSIS 

The objective of the analysis 1s to determine the radiation pattern of the 
antenna and its gain in the presence of a radome. In order to analyze the system. 
the spherical coordinates of igure 2.1 are used. The total electric field in the 
direction of unit vectors (8.0) Is required to compute the gain in the far field. The 
total electric field is the sum of the electric field due to the currents on the 
antenna (radiated field), and the electric field due to the currents on the radome 
(scattered field). The procedure is to assume a current distribution on the 
antenna, then use the MM to determine currents induced on the radome. In 
practice, scattering from the radome back toward the antenna will cause a change 
in the initial current distribution. This effect will be small for a well-designed 
radome and is ignored in the analysis. Scattering back toward the antenna will 
primarily affect the VSWR looking into its terminals. 

1. Determination of Gain 

Gain is a measure of the ability of an antenna to concentrate radiated 

power in a particular direction. The maximum value of gain for a lossless antenna 
is the directivity Go, which is the ratio of the maximum radiation intensity to the 


average radiation intensity [Ref. 3:p 610]. 


where Ula, 0 | is the radiation intensity and dQ = sin @d@do is the differential 


solid angle. Throughout this paper a lossless antenna will be assumed and the 





Figure 2.1 Coordinate System 


terms gain and directivity used interchangeably. 
Radiation intensity is the time-averaged power per unit solid angle, and 


IS proportional to the magnitude of the electric field squared 


u=s[f. 2-2) 
2n 


where 


E\= E*1+/E'|, (223) 














and E* is the electric field from the antenna, E° is the electric field scattered by 
the currents induced on the radome, and 7 is the impedance of free space. 
2. Determining the Electric Field 

The total electric field is calculated from the current distributions on the 
antenna and radome as illustrated in Figure 2.2. A current distribution J, on the 
antenna is assumed known. The currents induced on the radome are dependent 
on the current distribution on the antenna as shown in Figure 2.3, and this 1s 
taken into account in the method of moments solution. 

Maxwell's equations. are used to determine the electric field at a point in 
Space due to a current density J. When both the radome and antenna are 


present, the total current density 1s J=d,+d,. A solution of Maxwell's 


equations 1s 





E = -joA-——V[V oR) (2-4) 
where 
= = jkr 
A= 2 [fdlr.8.0) ds’, (2-5) 
and 
r=|R-R, (2-6) 


[In (2-4) Al is the vector potential, R_ the position vector to the observation point 


and RM’ the position vector to a source point on the antenna. 





Figure 2.2 Electric Field from Currents on Antenna and Radome 





E*(r,0,0) 


Figure 2.3 Electric Field Impinging on Radome 


B. ANALYTICAL TECHNIQUE 

Closed form mathematical solutions of the gain and radiation patterns are only 
possible for a limited number of applications that neglect the presence of the 
radome. In general, determining the radiation pattern and gain with the radome 
requires a numerical solutiton of Maxwell's equations. The radome geometry 
results in an unknown current distribution so (2-4) cannot be solved directly. 
However, boundary conditions can be applied on the radome surface which yield 
an integral equation for the current. The integral equation will be solved 
numerically using the method of moments. 

1. E-field Integral Equation and the Method of Moments 

The method of moments (MM) is a procedure used to solve an integral 

equation obtained by enforciny (2-4) on the radome surface [Ref. 5]. Equation (2- 
4) is used to determine the electric field £* incident upon the surface of the 
radome from the currents on the antenna J,. Initially the surface will be assumed 
a perfect electric conductor (PEC), in which case the total E tangent to the surface 
of the radome is equal to zero. Later a correction term will be added to the perfect 
conductor solution to account tor the dielectric properties of the radome. 


The boundary condition on the PEC radome 1s 


ee = Spe. (2-7) 
As applied to (2-4), 
Evan = Jo! {{ JG J y\y.(fd.ed (2-8) 
tan = JO || Gs + ° {| J.Gds |, 7 
; (WE ; 


where 9 


e~Jkr 





7 4nr 
and 


a R -R’). (2-6) 





Equation (2-8) is the E-field integral equation (EFIE). The only unknown is J. 
and once determined it is possible to calculate the scattered field by applying (2-4) 
to the radome. 


To solve (2-8), the current J, is expanded into a series, 


N = 
Jo = > Ib. (2-10) 


The J, are basis functions, and the I are the expansion coefficients to be 
determined. 
Mautz and Harrington applied this technique to the specific case of 


surfaces S generated by revolving a plane curve about the z-axis [Ref.5]. The 


plane curve is approximated by a series of straight-line segments connecting N, 
points. The surface has a circular cross section in the azimuth variable 0, thus the 


body is represented by a series of "faceted rings” as shown in Figure 2.4. 


10 


The surface coordinate system ts cylindrical [p, ¢,Z), and t is a vector 


tangent to the surface and orthogonal to 9 at every point. In order for the series 


expansion of the current (2-10) to approximate an arbitrary surface current density 


J,. independent sets of functions are defined as follows: 


Triangle Functions in t 





Pulse Functions in 0 — 


Figure 2.4 Segmented Radome with Current Basis Functions 
Depicted 








er eee , 
mt: J,, =t Ci a= Opa eee of 9, Cy cent None 
Ps 
= eg 
ino J, =0 Bs Hea ee See J=1,2,...,.Np-1. (2-11) 
Ps 


T, (t | is a triangle function extending over two segments, and P, it isa 


pulse function extending over a single segment. For an explanation of these 


functions see (Ref. l:pp 16-18]. Each value of f is referred to as an azimuthal 
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mode. The series expansion of (2-10) in terms of the chosen basis functions 


becomes, 
-3 By (a+ Sew (2-12) 
; ee a nj ng} ao 1 


In order for the series to be an exact solution, the sum must be extended over 
infinite modes -e<NM<oo. Of course the series must be truncated, and the 
maximum number of modes for convergence f,,,, 1S dependent upon the antenna 
scan angle and radome size. 


When the series expansion is applied to the integral equation (2-8), 


Es = J iw [Jou 6-=2(V dh, )06 as 


N=-2% sl 
e° =F 50 il Jou, 6 -=L(v" od, Jo is (2-13) 
PG |S WE 


In order to solve for the unknowns /,, and [,, > both sides of (2-13) are multiplied 


by testing functions Wi, and W° . and integrated. Using Galerkin's Method 


(Ref. 8], the testing functions are chosen to be the complex conjugates of the basis 


functions J’’’. 
ny 








Wi, = t au em 
Ps 

eet 

Ww’ =90 ew (2-14) 
Ps 
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Multiplying both sides of (2-13) by each of the testing functions and 


integrating results in a set of 2M, — 3 simultaneous linear equations for each value 


of m. In matrix form, these equations are written as 
[UV] = |Z]. (2-15) 


where U] is the excitation vector, [Z| IS an impedance matrix and 1 contains the 


unknown current coefficients. Submatrices associated with the t. and 9 


components (2-15) can be identified as follows: 


va | jh let Ge 
vey) (LZe"] [2a] be ~ 


The unknown coefficients are determined by solving the matrix equation 
a 
[W =(Z] [vy (2-17) 


The MM solution for the unknown current coefficients /,, and [? | in vector form 


is 


Z any (2-18) 


Once these are determined, it is possible to compute the electric field due to the 


current density distribution J, on the surface of the radome at any point in space. 


The elements of each submatrix of the excitation vector (U| are 
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Ue = ‘i WW? +E, ds, (2-19) 
and the elements of each of the submatrices of the impedance matrix [Z| are 


Zi | =|) el) Jous',6 - —= L(y 3, |v asas 


| 
Z|, = =i i, + |) Jn G a alt ree 
fn 


~ =1 Wj,» || ¥ = E (v’ ou’, \V 6 asas (2-20) 


The fact that the testing functions wi are orthogonal in 9 to Je due to 
the presence of the exponentials has been used to eliminate elements of (2-20) 
with n#mM [Ref. 6]. This allows each mode f to be treated independently, 
thereby simplifying the matrix inversion problem. Thus, by constraining the 
surface to be rotationally symmetric, the azimuthal dependence of the current ts 
expressed as a sum over exponentials (a Fourier series) and each mode fits 
independent. 


a. Impedance of Thin-Shell Dielectric Radomes 
The impedance matrix [Z| of the MM solution assumes that the 
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radome surtace is a PEC, where (2-7) is applied as a boundary condition of (2-4). 
For thin-shell dielectric radomes, a correction term is added to the impedance 
matrix that allows the modeling of radomes of arbitrary dielectric constant [Ref.9]. 

When the electromagnetic wave from the antenna is incident on a 


dielectric radome that has an intrinsic impedance 


oe W/ (2221) 
€ 


different from that of free space (all there will be reflection and refraction as 


shown in Figure 2.5 [Ref. 3:p. 397]. The total fields are 


satel 


— 


V x (H? + Hj = joe, (E +E*)+ jole - é, |(E° - ES |. (2-24) 


where the second term exists only when there is material present with a 


permittivity different than free space. This term is a polarization current J, where 


J = jole-e, JE. (2-25) 


ts 


Maxwell's equations are also satistied by the field from the antenna (E*, H? }. 
2 


Since (E, H| and ieee satisty Maxwell's equations, by superposition and (2-22) 


(E° , HS | also satisfy Maxwell's equations. Therefore 


V x HS = jok’ +d. (2-26) 





Figure 2.5 Discontinuity of Permittivity for Dielectric 
Radomes 


For a radome consisting of a thin dielectric shell, the polarization 
current has a large tangential component and small normal component. If only the 
tangential component is significant, then (2-8) can be modified to include 


dielectric properties of the radome as follows: 
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2d, 
Jjole - €, jf 


PR 


I] _ 
~ Jou] J,Gds+——V 8 [[a.eas (2-27) 


The first term on the right side ts the impedance loading term; the tntegro- 
differential portion remains unchanged from the PEC assumption. 


After multiplying both sides of (2-27) by the testing functions 
[Ws 


— WwW, and integrating. the 1mpedance matrix |Z| of (2-15) becomes 


[Z| oz (Aa + (Z, |. (2-28) 


where the elements of Zs ai are given by (2-20), and 


rae - | Wi oJS,,,Z,ds 

 ] 

Zu‘ |, = J LU! « oJ: ,2,d5 = 8. (by orthogonality of t. and 0) 
Zi), + | Ti . sleeeZ, ds = GQ, (by orthogonality of o. and t) 


Zr, - Jl We o3.,,Z,ds. (2-29) 


and 


a (2-30) 


Jjole ~ e, \t 


The thickness t in the impedance loading term is necessary to convert the volume 


current density J, of Figure 2.5 to a surface current density ds. 


7, 


Solving (2-17) yields 
1} = [Zan +Z,] [D]. (2-31) 


The radome impedance Z, can be rewritten in terms Of dielectric constanueaes 
thickness to wavelength ratio ,. and the loss tangent tan6 [Ref. 2:pp 45-46], 


7, aie 239 


} 


n,| dle, - 1} + e, tan 5] 


Note that both (2-29) and (2-32) allow the modeling of radomes whose 


permittivity may be complex, ¢ = €’ - Je” . 
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2. Mode Symmetry 
When the surface current density J, is expanded into a series as in (2- 


mee ie sui isvextended Over azimuthal modes N=@,+1,+2,--,+n... The 


mode index appears in the exponential term of basis functions J’ (2-11) and 


testing functions wr (2-14). The exponential factor can be written as, 
e/” =cos(no)+ Jsin(no). (23833) 


For each n there is a current vector {, given by (2-18). Mode symmetry refers to 
a relationship between the positive and negative mode vectors I, and _,. If mode 
symmetry exists, then only one matrix equation needs to be evaluated rather than 


two. This greatly reduces the computational effort. 


Recall the definition of the excitation vector [0] 


UP, =~ [ We, +E: ‘ds. (2-19) 


The electric field from the antenna 1s 


E? = £3(A,6,0)R+ E3(R,0,0)6 + E2(R, 0,0 )}0. (2-34) 


which on the radome surface can be expressed in terms of the two orthogonal 


components f and 
E? = E7 (t, o)f + E2(t, 09. (2-35) 
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In the far field of the antenna. the R component in (2-34) decays to zero and E? 


consists only of E%. Therefore as shown in Figure 2.6, 


E7 (t,o) = E7(t)coso 
E*|t,0) = E3(t)sino. 






Cut through sphere in 0= constant plane 


Radome Curve 


(2-36) 


Figure 2.6 Components of the Antenna Electric-Field on the Radome 


Plugging these far-field components into (2-19) 


*(t\p(tjat { cas o|cos(no) + jsin(no) do 


n 


oO 
w 
a 


ce 
li 
th ee, 
— 
=, 
a 
ss — 
Se 


~<— 
eS 


, 


2 es 


ll 
ee 
—_—| 
— —~ 
eee” | eee” 

oy 


=~ 
Il 
®@ 


o=6 
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(t)p(tjdt | sin o|cos(no) + jsin(no) do. 


(2-37) 


Only the © integration is dependent on nm. Expanding the integral in o and 


explicitly including the sign of gives 


cos o|Cos|+no | + jsin(+no }|do = 


2n 
o=8 
2X 
o=98 


2n 
COS 0 COS|+No|do + Jj | cos osin{+no|do. (2-38) 
o=0 . | 
where 
2° 
| cos osin(+no)do =i): (2-39) 
o=0 


due to the fact that the integrand is an odd function of @, and the integration is 


over one period () $ @ $ 2z. By trigonometric identity 


ii cos 9 COs(+No |do - ii cos|(+n- 1}o|do a ii cos|(=n+ Io|do. (2-40) 
0=8 =0 = 


o=8 


The right hand side is only non zero if M= +1. Furthermore, because cosine 1s an 


even function, the term U’, is independent of the sign of n. Thus 


Ue = ee (2th) 


For the second equation of (2-*7). the sign of U’. is dependent on the integral 


| sin o[cos(+n o}+ sin(+no) ido = 


o=@8 


MT) —.Y 


27 
" sino cos|=no jens in osin(+ no |do. (2-42) 


where 


big 
Pi deh o|do = 0. (2-43) 


by symmetry. By tmgonometric identity 


i sino sin{+no do = j cos|(+n z Io ldo = i cos|(+n+ Io|do. (2-44) 


0=0 o=8 0=8 


As can be seen in (2-44), the excitation vector U* is dependent on the sign of N 


due to the presence of the negative sign. Therefore 


| Tid | ae (2-45) 


In the near-field of the antenna the @ dependence is more complicated 
than the sine and cosine functions in (2-36). However, it has been verified 
numerically that ES and E? ave always even functions of @ and E% is always odd. 
This conclusion assumes that a symmetric amplitude and linear phase distribution 
excites the antenna. 


In light of the above results, the elements of the impedance matrix Zena 


of (2-18) have the following relationship between positive and negative modes: 


[Ref. 6] 


ho 
th 


Weg : tsa si us [Zemens |, |—Zusm ao (2-46) 


of 20 im of 90 
Zea |, Lawes) ~ Desatsay | Zeimron |, 


By inspection of (2-30), the elements of the load impedance matrix [Z| are 
independent of the mode N, because Wy and J, are complex conjugates of each 


other. Therefore. 


2) Zita . Hea UJ | — 


Rewriting (2-31) results in 


where 


I 


| = ae +Z,| . (2-49) 


By taking the inverse through partitioning [Ref. 8] it can be shown that the mode 


symmetry survives the inversion, and 


(an ey 


Liane ee 


By multiplying out (2-48), 


UUs, 


Equations (2-41), 


Die as Lellows: 


= i Do 


nj 
I" = = Weer y' + yor VU? 


mj A niyo nj? 


tf H 0 o ie t 
te lee 
oe at t 00 o — _fe 
av ty PN 0 ele 


(2-45), and (2-50) relate the positive and negative modes of (2- 


This is an extremely important result. Using mode symmetry reduces the number 


of calculations required to find the current coefficients ie by nearly half in most 


Cases. 


3. Measurement Matrices for the Scattered Field 
Once the current J, 1s obtained the scattered fietd E* can be determined 
by summing the contributions of the individual current elements in the far-field as 


shown in Figure 2.7. This is accomplished by use of a measurement vector with 


elements 
Rn = ||(R+6+0e” odds (2-53) 
FS 
In 
x far-field 
Z, 


Figure 2.7 Scattered Field 


In (2-53), a unit radiated plane wave is weighted by the current element at each 
point on the surface. Nezlecting the R component in the far-field, the 


measurement vector separated into @ and @ components is 


m 


Re = {J e Neer" Ba dds’ 
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R? = || geen —%™ 9 oJ, jds’. (2-54) 


5 


where U. U.and W are the direction cosines. 


U=SIN@COSO 
V=sinésino 


iv =Ccosé. (2-55) 


Thus, the components of the scattered field at a far-field point are the 


superposition of all surface segments 


— jk 1) 
ES = =k 1S Bint 


4rr 


pe alls Fin (2-56) 


4.  Closed-Form Far-Field Antenna Pattern 


When computing the antenna contribution to the total field in the far 


zone the higher oucletae terms can be neglected, and the radiation integral can be 


reduced to a closed form. The source of this field is assumed to be a uniformly- 


excited circular-aperture antenna scanned to a direction [8,5 05). Uniformly 


excited refers to a constant amplitude and a linear phase distribution of current on 


the antenna. The A, @ and @ components in the far-field are 


Fr =8 
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fe = =JRaCEa iv + nN, | 








47 rr 
_j ~jkr 
F? = & - (L, + nN,), Ca) 


where L,.L, N, and N, are given by 


[ = [[Mess gs’ = 0, (M= 8) 


, 


N = [[ dels“ ds’ = || (itd, + Gu, +Zd, Jer ds’ = (2-58) 


Ss’ 


where r’ and yw’ are as shown in Figure 2.8 [Ref. !O:pp. 455]. If the antenna 1s 


x-polarized, then N, = N, =@ and 


Ngee dom (2-59) 


Ss 


The equivalent in spherical coordinates 1s 


N, = || U, cosacosoe Ms" ds’ 


N, = || J, Sinoe*"°°5"'ds’, (2-60) 


s 


a 


The phase of the source point relative to the field point, r’COS w’, is a length that 


can be written in Cartesian coordinates as 


r’cosw =H'U+YU+2Z'W, (2-61) 


X 


Source Point l- 





Figure 2.8 Circular Aperture Antenna Source Point 


where U,V, and W are the direction cosines (2-55). Including beam scanning 


adds an exponential phase factor to the current 


7 Jk(U,+YU,+ZW, | 
ee 
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where U,. U, and W, are the direction cosines of the scan angles [8,5 0, which 


also obey (2-55). J, is the magnitude distribution, assumed to be constant, and 


OY Z’) are the coordinates of a source point. For an antenna iying in the x-y 


plane, Z’=@. Therefore, N, can be rewritten as 


k| U,+ U; —- ki s + s P 
Nees COs COSio rer ek Maumee 
E 


Jk(K'\u,-uj\+y'(v 


= J, cos#cosolle Ng? (2-63) 


By setting #’ = p’COS@’ and Y’ = p’ Sino’, the exponent term of (2-63) become : 


kp'|cos o’(u, -u)+sino’(v, -v)], (2-64) 


Defining the variables 


and considering these to be orthogonal vectors as in Figure 2.9, it is possible to 


rewrite (2-64) as 


oA ATSB" | COs ee sini ee 
VA2 + B VA? +B? 


ya 


From Figure 2.9 the factors containing A and B in the brackets can be replaced 


with COS y and siny 


Kp’\ A? + B’[coso’ cosy +sino’sin 7}. (2-67) 


By trigonometric identity. 





Kj)’y AC + B*|cos|o’ ~ 7) (2-68) 





Figure 2.9 Representation of A and B as vectors 


Using this form of the exponent in the integral of (2-63), gives 


2 td ae a 7. 
N, =u, Cos @coso | | Coa Bt c05(°'— 1} 4 do’ do’. (2-69) 


0’=8 p’=8 


which conforms to the definition of a Bessel function. The same integration 1s 


performed tn the second equation ot (2-60). and the same Bessel function results 


Finally, 


J\{ kay AP +B? | 
N, = J, COS O8COSo za*| 2 —__—~ 


ka\A* + B° 


| Kav A? + B? | 
N=. SINC ca ee (2-70) 
° ‘ kav A2 +B? 


where d, 1s a Bessel function of the first order. The far-field equations (2-57) 


become 
|) 2u,( kava? + B? | 
Jkr 1 
Poa= =U UNigeme =s7cos >. 40) —— 
ate kav A’ +B? 
kn em || 2u,( Kava? +B? | 
P, =e So... (2-71) 
4n or kavA? + B? 
where 
A=u,-u=siné, coso, -singcaso 
B=u,-U=sind, sing, -SINGsiNo, (are) 
and 


VA +B? = 





sin’ g+sin’ 6, -2sin@,sin@cos(o-o0,). (2-73) 


Using equations (2-71) through (2-73) are computationally less 
intensive than a numerical integration. These equations were employed for both 
the far field pattern and gain computations. 

5. Near Field of the Antenna 

In the previous section, the far-field of a circular aperture antenna was 
determined. When a radome is located in the near-field, the integrals cannot be 
reduced to closed form and they must be evaluated numerically. The near field is 
generally considered to be the region surrounding the antenna where Ff << A, or 
kr << 1([Ref. 10:pp. 106]. 

The complete expressions for the near-field of a circular aperture 


antenna excited by a current polarized in the x-direction are 





ies pe a x’) God, le away 








Ey =f |J||4-# )[Y- 9), PM aay 
E, = ae [lle ee a a 
where 
k‘r? -1- jk 
;- 
6, = 5 SJR (2-75) 


and 


— 


2 


2 pa | 
Pao — Hae Oe ie 2a | (2-76) 
The primed coordinates designate source points and the unprimed an observation 


point on the radome as in Figure 2.10. (Ref. 11] S’ is the area of the antenna 


aperture. Since the radome 1s rotationally symmetric, the components (one 


can be described in terms of cylindrical coordinates using the transformations 







[H’, y’,@) or ier Oo’, 8) 






P(x,Y,Z] or P(R, 6, 0] 


Figure 2.10 Rectangular Coordinates for Circular Aperture Antenna 


HK’ = p’COso’ 
y’ =p’singd’ 
dx’dy’ = p’dp’do’. (2a) 


For an antenna scanned to (8,0, | the current distribution 1s 


jk| #’coSs\ 0’ -0 sine +Y’sin| 0’ -o sing | 
J, =d,(p’ Je 2 . 2 >. (2-78) 


Applying (2-75), 2-76) and (2-77) to equations (2-74) results in 


2r a 
=; ae Kp’COS(o-0,|sine@, _ | 
le jn | fa | | Jkp’COSs|o-o0,}sin jkr 7 


0’=8 p' =8 


6 +(pcoso- p’coso’) G 2 pd ao 





a =i i fu, (0 ‘Je el cos o- os | ea jkr . 


o’=8 p’-8 


[| pcos 0-p’coso’)(psing -p’sin 0’), |p’ dp’ do’ 





ox a 
_ <i 1\  Jkp’ cos(o-o, )sing, 4- jkr 
aps = rn J Jel ie Bae AS aa 
(p COSo-—p’cos 0')Z6, |p’ dp’ do’ | 
where 
a (p-p’) +(2-2'} +4pp’sin?(*), (2-80) 


Since the expression for the excitation vector requires spherical components, 


(enn l2.| are converted using a coordinate transformation 
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F, =£, singcoso+E£, singsino+&, casé 


£, COSPCOSe+E, COSOsing—E, sing 


0] 


(2-81) 


E 
io Sino + EA COSo: 


a) 

These are the equations programmed to compute the E-field impinging 
on the radome. The magnitude of the current distribution is allowed to vary as a 
function of p’. Various antenna sidelobe levels can be modeled by simply 


; 


changing the function J, ("| 


Ili. PROGRAMS 


The mathematical formulas developed in the previous chapters were 
computer programmed using FORTRAN. The pattern and gain calculations are 
done in two separate programs: RADOME.F and GAIN.F, both of which run on a 
Sun SPARCstation™ under UNIX. A brief description of each follows. 


A. DESCRIPTION OF THE PROGRAM RADOME.F 

This program, developed by Francis [Ref. 2], determines the series expansion 
coefficients of the current J, induced on the surface of the radome. It also 
calculates the total electric field surrounding the system, and generates radiation 
patterns using MATLAB™. Originally called LDBOR.F, (an abbreviation of 
“Loaded Bodies of Revolution”) in Ref. 2, it has been modified to take advantage 
of symmetry, and also allows the antenna mainbeam to scan in both @ and o. 

Figure 3.1 isa flow chart of the program. To run the program, data 1s entered 
interactively, except for a data file “gaus/gaus(ITH)" containing the integration 
constants for Gauss Quadrature. which is read in the program. All integrals are 
evaluated using this technique (see Appendix A.). Table [ is a sample input. 

After all data is entered, the coordinate points of the surface are computed 
and stored in the arrays RH(= p) and ZH(= Z). The surface impedance of each 
segment is stored in ZLO. The elements of the MM impedance matrix are 
computed in subroutine ZMAT, and the excitation vector elements in GENEX. 
These provide par and U| in equation (2-17). Subroutines DECOMP and 
SOLVE solve the matrix equation using Gaussian elimination. Subroutine 
PLANE calculates the measurement matrices, which are subsequently used to 


compute the scattered field from the radome. 
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Two antenna radiation subroutines are included. CIRCRTP calculates the 
electric field components &,.£,. and E, according to equations (2-81). It is used 
by subroutine GENEX to obtain the excitation vector. The second antenna 
subroutine ANTFF assumes that the observation point is in the far-field. This is 


used in the main program to compute the antenna contribution to the radiation 


— 


pattern E® in (2-3). 

Output includes an ASCII file “outldbor" which lists all the input parameters, 
radome geometry and impedance. and the pattern magnitude. The current 
coefficients are written on a file “curcoefsdat” so that further patterns can be 
generated without recomputing the coefficients. Finally, MATLAB™ formatted 


files are generated for plotting the @and @ components of the electric field. 
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TABLE I. SAMIPLE INPUT RADOME.F 


ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS 


SELECT BOR GEOMETRY BY NUMBER 


NUMBER BOR 
I OGIVE 
2 ialleleds 
5, CONE 
+ DISK 
5 PARABOLA 


ENTER SURFACE CURVATURE (wavelengths) 


ENTER ZPRIME,WHERE CURVATURE STARTS (wavelengths) 


ENTER BASE RADIUS (wavelengths) 

ENTER FILENAMES gaus### 

ENTER THE FILENAME IN 7 (NT) 

ENTER THE FILENAME IN PHI (NPHD) 

ENTER THE FILENAME IN X AND Y 

ENTER THE PLOTTING INCREMENT IN DEGREES 
ENTER THE HIGHEST MODE 

ENTER PHI (observation) IN DEGREES 

ENTER SCAN ANGLES IN DEGREES: THETA, PHI 
ENTER COMPLEX IMPEDANCE, OHMS: (Real, Imag) 
ENTER THE ANTENNA RADIUS (wavelengths) 
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A 


iS 


90 

30,90 
(O,-1 700) 
l 


Begin 





Assign Parameters 


r. [ouput written to files appended by ITH 


Enter Geometry Selection=? 





1=Ogive J=oplere 3=Cone 4=Disk S=Parabola 


Y 


Call PARAB 
Determine NP. 
ZH,.RH,.DM, 

FOD,ZP 


Call OGIVE Call CONE 
Determine NP, Determine NP, 
Zi RH BASE. ZH.RH.BASE. 
RS,ZP RS.ZP 






Call SPEIERE 
Determine NP, 
ZH.RH.BASE, 
RS.ZP 


Call DISK 
Determine NP, 


ZR OAS. 
Ro Ze 












NP= # poimts on curve generating BOR 
ZH,RH=Coordinates of BOR 
DM=Diameter of Parabolic Radome 
FOD=r/> 

ZP=Location of antenna z' 
BASE=Base diameter of antenna 


Enter filenames for Gauss in t 
and © 


Figure 3.1 RADOME.F Flow Diagram 





Read XT,XP,X,A.AT,AP 


a | 
<a‘ 
Enter DT, Plotting Increment 
Yes 


Calculate NBLOCK.NROW 


NBLOC partitions current vector [I] 
into by mode n 

NROWs=total # of individual elements 

of [I] 





Weights and Abscissas are read in from 
files named Gaus/Gaus(ITH) 




















\ 


Enter P, Observation Angle 
No <p 
es 


Enter THD,PHD,IMP,ARAD 


Read ITH,SELECTION.BASE 
NP.MT,.MP,N,NROWS, 
MODES.NBLOCK,.THS.PHS. 
ARAD,PHI,IMP,ZP.RS,RH. 
ZH.ZLO from "curcoetsdau TTH)" 









Calculate IT.MHILMP,MT,N 
= = Prepares tor E-field calculation 


Figure 3.1 RADOME.F Flow Diagram ( Continued) 
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Exit 







: Loop 
pes2— NP 


Calculate ZHB,.RHB.ZLO() 


Assign surface impedance to each 
segment of radome 













Yes 
Call ZLOAD | 
rewlinseZ5 
Subroutine ZLOAD calculates 
= elements of (Z; ] 
Exit Loop 


DO 400 MelMHI = 


; -_ — = Positive mode n 


: Call ZMAT 
etums Z 
eee Subroutine ZMAT calculates 
elements of (Zany) 
Call ZTOT | 
retums Z 
z Subroutine ZTOT calculates 
7 elements of (Z]=(Z yy }+(Z_ | 
Call GENEX 
retums B 









Subroutine GENEX calculates 
elements of excitaion vector [V ] 


Call | 
‘et Z 
ee Subroutine DECOMP calculates 


Figure 3.1 RADOME.F Flow Diagram ( Continued) 





elements of (Z,y! 


4] 


Call SOLVE 
retumms C : 
Subroutine SOLVE calculates 
elements of [IJ=[Z] -![V] 
Store CT, (+n) current 
elements 


— = Negative mode -n 





Call GENEX 
Call ZMAT 
Call ZTOT 
Call DECOMP 
Call SOLVE 








Write | THSELECTION BASE,.NP.MT,.MP,N 
NROW .MODES,NBLOCK,THS,PHS,ARAD. 
PHI IMP.ZP.RS .IFF.RH(L).ZH(L).ZL(O),ZLO(), 
GNT.C:NPHI.CGP to curcoefsdat(ITH) 


These parameters are to be used 
“a in Gain computation 





Figure 3.1 RADOME.F Flow Diagram ( Continued) 
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Exit 





Compute E-field for each angle 
ot plotung increment 


Loop 
DO =i 


Initialize ET!.ET2 
ErPiBP2 


Exit Loop 
DO 300 M=1.MHI 


Calculate EXP1.EXP2 
: et/-jnd (2-56) 


Call PLANE 
retums R 








Compute scattered eld) E*(2-13) 















Subroutine PLANE calculates 
elements of Measurement Vector 
R (2-54) 






Loop 


DO 402 L=lMyF 








Calculate ET1,EP1 


ETI.EP1 are elements 





Loop 


DO 403 L=MT+1,N 
Calculate ET2,EP2 


ET2.EP2 are elements 






< 


Figure 3.1 RADOME.F Flow Diagram ( Continued) 









Call ANTEP Call CIRCRTP returns 
returns ETF, EPF CIRCR, CIRC, CIRCP 


Subroutine CIRCRTP 
claculates near field of 
antenna 








Subroutine 
ANTFI- 
calculates 
far-field 

of antenna 










Calculate ETF.EPF 
from CIRCR, CIRCT,CIRCP 













Calculatc ER RESEPEE. 6 fee 
ETSGAT. ERSGATAEC EX 
BC VEX Vee hale \R: 
EXI, Pii@ePHxXyANG: EC x 






(8:0) components of E-field, total field. 
scattered field, copolarized, cross polarized, 
phase, angle 


Write E-field components to 

separate files 
E-field files are formatted for input to 
MATLAB | for plotting purposes 


Figure 3.1 RADOME.F Flow Diagram ( Continued) 
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B. GAIN.F 

This program determines the gain of the system by solving (2-1). In order to 
pertorm the integration, the current coefficients from RADOME.F must be 
provided so that the radome scattered field can be computed. Input consists of 
specifying an iteration letter, and selecting the number of divisions in the @ and o 
integration. Input data is then read from a file “curcoefsdat" that is appended 
with the iteration letter chosen. 

Figure 3.2 is a flow chart of the program, and Table II is a sample input. 
GAIN.F contains several of the same subroutines that occur in RADOME.F. They 


include PLANE, CIRCRTP, and ANTFF. 


TABLE Il. SAMPLE INPUT GAIN.F 


ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS 
ENTER NUMBER OF DIVISIONS IN THETA DIRECTION, 


NeIvi =? 
ENTER NUMBER OF DIVISIONS IN PHI DIRECTION, NDIVP =? 
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Assign Parameters 















Read [TH.S]-LECTION, BASE 
NP.MT.MP.N.NROWS, 
MODES.NBLOCK.THS.PHS, 
ARAD,PHI.IMP.ZP.RS.RH, 
ZH.ZLO trom “curcoefsdat(ITTH)" 


data was written to "curcoefsdat(ITH)" 
from RADOME.F 


CT(L) are elements of current coefficient 
vector [I] written to curcoefsdat 


Readewdid), ATU). ACD ex P( CAR 


A(T) are weights, and X(I) are adscissas 
which are read from gaus/gaus# 



















NDIVT, and NDIVP are the number of 
divisions in @ and @ respectively. This 
allows for easy variation of the number 

of integration steps without recomputing 
weights and abscissas for Gauss 
integration. 








Calculate S(1). 6 integration angles 
and Q(1), ® integration angles 


Set SUM, SUMNONE, UMAX, UMAXNONE, 
EMAX, EMAXNONE to 0 





Figure 3.2 GAIN.F Flow Diagram 
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Exit Loop 
Do 300 JJ=1.NDIVP = Summation of divisions in © 






Do 300 Jai NPH! 


Summation within each division ino 







Loop 


Do 301 1l=1,.NDIVT 


Summation of divisions in @ 









Summation 
within each 
division in @ 


Calculate THR.PHR 


Loop 


Do 30! lai,NPH! = 













E-field determined at 
each THR (8), and 
PHR (®) angle 


Seee eer ret wall, EP2 to.0 


Exit 


i 
Do 305 M=1.MHI se 






Mode, n. Summation 


of E-field (2-13) 
Calciilatereil 242 
et/-jngd (2-56) — ~ 


Call PLANE 
(6) returns R 
Subroutine PLANE calculates 
C+) >) & elements of Measurement Vector 
R (2-54) 


Figure 3.2 GAIN.F Flow Diagram ( Continued) 
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©) 


DO 306 L=1.MT 
Calculate F T.EP Calculate ET1,EP1 


ETI,EPI are elements 
of ES, (2-56) 














Calculate 
Eiiser 


xa Loop 


Xit 
DO 307 L=MT+1.N 
Calculate ET2,.EP2 
<> No 


ET2,EP2 are elements 






of ES, (2-56) 
Call ANTFEF Call CIRCRTP retums 
returns ETF, }:PF CIRCR, CIRC, CIRCP 
Subroutine CIRCRTP 
claculates near field of 







Subroutine 
ANTFF 
calculates 
far-field 

of antenna 





antenna 
Calculate ETF,EPF 
trom CIRCR, CIRCT,CIRCP 


Calculate ETCOMB, EPCOMB 
EMAG, EMAGNONE. UMAG., 
UMAGNONE, SUM. SUMNONE Total E-field in 8 and 0, magnitude ot E-field 
with and without radome, magnitude of field 
intensity with and without radome, current 


sum with and without radome. 






Figure 3.2 GAIN.F Flow Diagram ( Continued) 
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Calculate PRAD, PRADNONE 








PRAD= radiated power with 
radome 

PRADNONE= radiated power 

without radome 






Calculate G.AIN,. GAINNONE, 
GIB. GDBNONE 





GAIN= gain with radome 
GAINNONE= gain without 
radome 
GDB=gain with radome, (dB) 
GDBNONE=gain without 
radome, (dB) 





Write gain values to “gainout(ITH).m" 


"gainou(ITH).m" is formatted for 


input to MATLAB !M 
for plotting purposes 





Figure 3.2 GA/N.F Flow Diagram ( Continued) 
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IV. RESULTS AND CONCLUSIONS 


A. VALIDATION AND CALCULATED RESULTS 

Part of the continued development of this research 1s validation of the 
programs. Test cases were run to see how the programs performed against 
known results. Test radomes ranged from nearly nonexistent disks located at 
large distances, to perfectly conducting parabolas located in the near-field. ie 


plane and H-plane patterns defined in Figure 4.1 were plotted to determine the 


E-plane H-plane 


@ = @° Plane 


TKK bg Uo 
LZ 


Figure 4.1 Orientation of E-plane and H-plane 





. 


~, 


=~ 


A 
VL” 


etfects of the radome. The antenna pattern without the radome was plotted 
(dashed line) along with the pattern of the complete system (solid line) to 


illustrate the effects of radome scattering. 


50 


1, Code Validation: Scattering from a Small Disk 
The pattern of an antenna radiating in the presence of a negligibly -small 
circular-disk located at a large distance was computed. Since the radome 
scattered field is negligible. the predicted output is that of a uniformly-excited 
circular-aperture antenna. Figure 4.2 shows the plots of the E-plane pattern for 
three different radius antennas (.5 4, 1.04 and 2.02). The plots appear as single 
lines with no scattering. The gain as computed by GA/N.F. was compared to the 


theoretical gain of a circular aperture with a uniform current distribution 


(4-1) 








Where @ is the aperture radius [Ref. 10:pp. 483]. Theoretical gains for the .5 A. 
1.0 A and 2.0 A radius apertures are 9.94. 15.96. and 21.98 dB respectively. [n 
each case the numerical gain is within |.5 dB of the theoretical gain. 
2. Code Validation: Perfectly Conducting Paraboloid 

Programs have been written that use the MM solution for bodies of 
revolution developed by Mautz and Harrington to predict the patterns and gains 
of reflector antennas [Ref. 7]. Assuming that the radome ts a perfects 
conducting parabola. the output of RADOME.F and GAIN.F can be compared 
with the output of these existing programs. 

Figure 4.3 is a test case where the radome is a parabola with surface 
impedance Z, = @+ J@ Q. The dimensions of the radome matched those of a 
reflector which was analyzed by a separate computer program that calculates 


both pattern. and gain. The pattern and gain for the perfectly conducting radome 


oy 










Antenna Radius = .5 


Numerical Gain = 11.42 (dB) 
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Figure 4.2 Scattering from a Negligibly Small Disk 
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computed by RADOME.F and GAIN.F are identical to that computed for the 


reflector. 


Perfectly Conducting Paraboloid 








Diameter=5A 
Focal Length =2A 


Feed 






Backscatter from Reflector —O 
be 


Gain = 23.04 (dB) 





50 100 150 
Theta Angle (deg) 


Figure 4.3 Perfectly Conducting Paraboloid 


3. Dielectric Radomes 
To determine the effect ot various radome materials on antenna 


performance, radomes with several different surface impedances were analyzed. 


Surface impedances were calculated using (2-32) for a range of e€, and tané 
typical of radome materials. The real and imaginary components of Z, as a 
function of €, and tano are plotted in Figure 4.4. As shown in the figure, the 
reactance 1s essentially independent of tand over the range O<tan6d <@.1, 
and the real component of Z, 1s strongly dependent on tan, particularly in the 
region of low é,. 

A test radome consisting of an ogive with dimensions given in Table I 
and parameters 1, = 1/28 and tan6 = 0.01 was analyzed.* The patterns were 
calculated for dielectric constants ranging 2<e, <1@, and the results are 
summarized in Table III. The E and H-plane patterns of an ogive radome with 
Z,=34- Jj1700 © are shown in Figures 4.5 and 4.6. The impedance 
corresponds to a radome made out of a material of ¢, = 2.0 and tan6é =@.81. 
The scattering 1s symmetric with respect to 9 for the unscanned antenna, and the 
backlobe level is -45 dB below the peak. The scanned case has higher sidelobes, 
as expected, with a lobe of -30 dB at an angle 270°. The gain loss due to the 
radome for the unscanned case 1s approximately 0.1 dB, and for the scanned case 


0.2 dB. The discontinuity in the H-plane patterns at @ = 9@° is due to the 


antenna model. When @ = 9@° the field is @ polarized, and according to (2-71), 
E* has no COS@ factor to force the pattern to 0 at @=9@°. Note that the 
radome ohmic loss has not been included. 


Increasing the dielectric constant results in a constant-shape scattering 


pattern with increasing magnitude. Figure 4.7 is the E-field pattern for e, = 10, 


which corresponds to some ceramics. It is apparent that the scattering by the 


" These dimensions are not intended to match any particular existing 
radome. 
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radome in the rear hemisphere ts significant, and the gain loss is nearly 4 dB. 


Table If] shows the gains with and without the radome for 2 < e, < 10. Figure 
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Figure 4.4 Surface Impedance Z, as a Function of ¢, and tano 


4.8 is a plot of the gain loss relative to no radome as a iunction of é€,. For 


reference the Fresnel reflection coefficient for a planar interface 1s also shown. 


TABLE III. GAIN LOSS FOR AN OGIVE RADOME 


Z, (2) 


ual 
28” 


tand =@.01 


34-1700 
7.5-7560 
5.3-j420 
3.3-j280 
2.7-j240 
2.4-j210 
2.1-j190 


| 


Gain (dB) 


With/W ithout 


Radome 


6, = 0° (0, = 8°) 


Phe 


e2e3 
22s 
BS, 


21.42/22 
20.91/22.32 
19.95/22.32 
19.51/22.32 
19.04/22.32 
18.59/22,32 


2 


5:6 


Gain (dB) 


With/Without 


Radome 


@, = 50° | o eae 


ZILGI/2 Ts 
2-6 1/2 is 
20,90) lees 
20.46/21.78 
19.47/21.78 
18.96/21.78 
18.47/21.78 
18.00/21.78 
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Figure 4.5  E-Plane Pattern of Antenna with an Ogive Radome 
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Figure 4.6 H-Plane Pattern of Antenna with an Ogive Radome 
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Figure 4.7. E-Plane Pattern of Antenna with an Ogive Radome 
Z,=2.1- jf187 Q 


7) 
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The gain loss appears to be linear with €,. This is probably because the base of 
the ogive 1s open, which allows reflected waves to exit the radome. If a back 


surface was present, these waves would be reflected again causing additional 


gain loss. 
a. Conical Radome 


To determine the effect that radome shape has on scattering, a 


conical radome was analyzed. The length, width, and impedance of the cone was 
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Relative Dielectric of Radome 


Figure 4.8 Gainloss Relative for an Ogive Radome as a Function of ¢, (E- 
Plane Scan) 


the same as the ogive radome discussed previously. A conical radome presents a 
relatively flat surface that could possibly cause more focused scattering than a 
curved shape like an ogive. Ray tracing of Figure 4.9 shows that the reflections 


from the radome will add coherently for the cone, and non-coherently for the 
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ogive. However, this assumes that the radome is in the far field of the antenna. 


The E-plane pattern is shown in Figure 4.10. 


Main Beam Refraction 





Scattered Field 





Figure 4.9 Scattering from Conical Radome vs. Ogive 


A comparison of Figures 4.5 and 4.10 demonstrates that there 1s 
even less scattering from the conical radome in the rear hemisphere. For the 
scanned case, the conical radome has a distinct lobe at 6 = 28@° due to the 
specular reflection of the main lobe from one side of the radome. This is referred 


to as a reflection lobe. 
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Figure 4.10 E-Pane Pattern of Antenna with a Conical 
RadomeZ, = 34 - jl 788 Q 





B. CONCLUSIONS 

This research has resulted in two computer codes that can be used to predict 
the pattern degradation and gain loss due to an axially symmetric radome in the 
near field of a linearly polarized circular aperture. The computed patterns show 
that the gain loss is slightly greater than the reflection loss for a planar interface 
with the same dielectric constant as the radome. This is approximately true for 
low-loss radomes. For lossy radomes the impedance has a very strong affect on 
scattering, while radome shape has less effect. Scanning the antenna results in an 
asymmetric scattering pattern and generally higher sidelobe levels. 

The computed gain compares favorably with values computed by closed form 
solutions derived from theory, and the gain of the system has the expected 
decrease as the level of scattering increases. For the specific case of the ogive 
radome, there is a linear dependence of the gain loss to the relative dielectric of 
the radome. Note that the gain calculations presented here do not include ohmic 
loss inside of the radome. 

Improvements in the original code include arbitrary aperture illumination and 
incorporating mode symmetry. Computational run time is reduced by nearly halt 
when symmetry is used, and significantly increases the applicability of this 
method of analysis. 

In order to apply this analysis to multilayer radomes, an equivalent 
surface impedance must be determined. The program can be modified to account 
for radomes constructed of nonuniform materials by assigning different 


impedances for each segment of the plane curve generating the BOR. 


63 


APPENDIX A. GAUSS QUADRATURE 
Computing the MM elements of Uj and [Z| (2-19 and 2-28) as well as the 


gain requires the numerical evaluation of integrals. MM requires numerous 
integration in order to solve for the current coefficients /,, and f, of (2-3 er 
only is numerical analysis necessary to perform the integrals in the MM. but also 
the matrix analysis of (2-31). The large number of separate integrals necessitates 
an integration technique which provides accuracy in a few steps to avoid 
excessive run time. Gauss Quadrature is a method of numerical integration which 
employs unequally spaced intervals , and provides accuracy with a few points. 
[Ref. 4:pp. 154-157] 


Gauss Quadrature approximation of the form 


b 
l= | f(x)dx., (A-1) 
by the sum of m terms 
P= 225 wf, ). (A-2) 
ce —) 


The coefficients W, are weights, and the #,'s are abscissas determined from 


b+a b-a 
= ate : (A-3) 
k 5 5 Cig 








where &,'s are the m number of zeros of the m'th degree Legendre Polynomials. 


APPENDIX B. RADOME.F LISTING 


The material on the following pages is a listing of the program RADOME.F 


described in chapter IIL. 
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CC a a RR KO OO KE 


C ROGRAM > radome.f Views 


C DATE = 25 Jantiary 1992 
C LATEST REVISION: 31 December 1992 
C PROGRAMMERS : D. JENN, R. FRANCIS, K. KLOPP 


C 

C NOTES: 

C 1. ARBITRARY CIRCULAR APERTURE ILLUMINATION SPECIFIED IN 
C FUNCTION TAPER. 

C NEAR FIELD COMPONENTS INCLUDED. (COMPUTED IN SUBROUTINE GENEX. ) 
C ANTENNA IS AT THE BOR COORDINATE SYSTEM ORIGIN. 

C LINEARLY (X-POLARIZED) ANTENNA. 

C CLOSED FORM FAR FIELD ANTENNA OPTION ADDED. 

C FLAGS USED THROUGHOUT: (IMP IS SET AUTOMATICALLY) 

C IMP=0 perfect conductor radome (for test purposes) 
C IPRINT=0 print pattern points to unit 8 

C ISEG=0 print the generating curve points to unit 8 
C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


On fF WN 


ICALC=0 compute current coefficients & pattern 
#0 read current coefficients from disc file given 
by ’CURCO and compute the pattern. 
IFF=0 closed form far field antenna pattern is used to 
compute ETF and EPF (SUBROUTINE ANTFF) 
#0 far field computed using CIRCRTP at distance ROBS 
ISYM=0 use mode symmetry 
#0 turn off mode symmetry (for testing) 
7. THE FACTOR ETA (=377 OHMS) IS OMITTED THROUGHOUT. 


BEGIN MAIN PROGRAM¥## KEE EEE KEKE KEKE EERE KEKE KEKE EERE 
CHARACTER ITH 
CHARACTER*9 OUT,ETSCATTER,EPSCATTER,  CURCO 
CHARACTER*12 CC 
CHARACTER*6 AN,ETFEED, EPFEED 
CHARACTER*8 GNT,GNPHI,CGP,CPOL, XPOL 
CHARACTER#14 TPTS,PPTS,PHIPTS 
COMPLEX EP,ET,Z(100000) ,R( 1600) ,B(800) ,C(800) ,U 
COMPLEX UC,ET1,EP1,ET2,EP2,EC,EX,ZL0(400) ,ZL(2400) 
COMPLEX CIRCR,CIRCT,CIRCP ,ETF ,EPF 
COMPLEX EXP1,EXP2,CONJG,CEXP,CMPLX ,CT(30000) , IMP, JK 
DIMENSION RH(400) ,ZH(400) ,XT(4) ,AT(4) , IPS(800) ,XP(100) , AP(100) 
DIMENSION A(100),X(100),EXP(500) ,ANG(500) , ECP(500) 
DIMENSION ECV(500) ,EXV(500) ,PHC(500) , PHX(500) 
DIMENSION ETSCAT(500) , EPSCAT(500), ETHF(500) , EPHF(500) 
INTEGER CNPHI (SELECTION 
DATA IPRINT/O/,ICALC/1/,IFF/0/,ISEG/0/, ISYM/0/ 
DATA START,STOP/O. ,360./,PI/3.1415926/ 
Rad=PI/180. 
ECK=0. 
BK=2.+*PI 
U=GOre i) 
UO=(0.,0.) 
UC=-U/4./PI 


JK=(0.0,6.283185307) 
C*xx«x«e*e ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS *«*xex «naxx 
C#* ex (ALL DATA FILES ARE APPENDED WITH THIS LETTER) KR KK KOK KK 
Cee Ke (.m FILES ARE FOR GENERATING PLOTS IN MATLAB) KKK RK KH 
WRITE(6,*)’ ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS’ 
READ(S,*) ITH 
OUT=’outldbor’//ITEH 
AN =’ang’//ITH//’ .m’ 
BEPEED=*ett’//118//2 :m’ 
BEPEED="epi 7/711H//* .m’ 
ERSCALIER= seteceat//L1H//7?.m 
EPSCATTER=’epscat’//ITH//’.m’ 
CURGO="asciidat //11H 
CG="eurcoersdat’ //11H 
CPOL='etpol’//ITH//’ .m’ 
XPOL= eppol’//1TH//* .m’ 
IFCICALC.EQ.0) THEN 
WRITE(6,*)’ENTER THE ANTENNA RADIUS (wavelengths) ’ 
READ(S,*) ARAD 
C*****eSELECT THE GEOMETRY OF THE BOR##X& Ke 4K KKK EEK REE 
WRITE(6,*)’SELECT BOR GEOMETRY BY NUMBER.’ 


WRITE(6,*)’ NUMBER BOR’ 
WRITE(6,*)?’ 1 OGIVE’ 
WRITE(6,*)? 2 SPHERE’ 
WRITE(6,*)? 3 CONE’ 
WRITE(6,*)’ 4 DISK’ 
WRITE(6,*)’ 5 PARABOLA’ 


READ(S,*)SELECTION 
IF (SELECTION.EQ.1)THEN 

WRITE(6,*) CALLING SUBROUTINE FOR OGIVE’ 

CALL OGIVE(NP,ZH,RH,BASE,RS,ZP) 
ELSEIF (SELECTION .EQ.2) THEN 

WRITE(6,*) CALLING SUBROUTINE FOR SPHERE’ 

CALL TESTSPHERE(NP,ZH,RH,BASE,RS,ZP) 
ELSEIF (SELECTION .EQ.3)THEN 

WRITE(6,*)’CALLING SUBROUTINE FOR CONE’ 

CALL CONE(NP,ZH,RH,BASE,RS,ZP) 
ELSEIF(SELECTION.EQ.4)THEN 

WRITE(6,*) ’CALLING SUBROUTINE FOR DISK’ 

CALL DISK(NP,ZH,RH,BASE,RS, ZP) 
ELSEIF (SELECTION .EQ.5) THEN 

WRITE(6,*) CALLING SUBROUTINE FOR PARABOLA’ 

CALL PARAB(NP,ZH,RH,DM,FOD,ZP) 

BASE=DM/2. 
ENDIF 
IF (NP.GT.400)THEN 

WRITE(6,*)’ MAXIMUM NUMBER OF POINTS(NP) IS 400’ 

GOTO 999 
ENDIF 

ENDIF 


CSP ESE SESS SESE SSS ESSE TS LES ESTES ESSE SSE SESE SES SESS EES SE SESE SS 
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WRITE(6,*)’ ENTER THE FILENAMES gaus###’ 
WRITE (C6 Jes poco eee , 
WRITE(6,*)’ ENTER THE FILENAME IN T (NT)’ 
READ(5,*)GNT 
WRITE(6,*)’ENTER THE FILENAME IN PHI (NPHI)’ 
READ(S , *) GNPHI 
WRITE(6,*)’ENTER THE FILENAME IN X AND Y’ 
READ(S,*)CGP 
C OPEN THE FILES FOR THE gaus/gaus### 
TPTS=’gaus/’//GNT 
PPTS=’ gaus/’//GNPHI 
PHIPTS=’gaus/’//CGP 
OPEN(1,FILE=TPTS ,STATUS=’OLD’ ) 
OPEN(2,FILE=PPTS ,STATUS=’OLD’ ) 
OPEN(4,FILE=PHIPTS ,STATUS=’OLD’ ) 
READ(1,*) NT 
IF(NT.GT.4) THEN 
WRITE(6,*)’MAXIMUM NUMBER OF POINTS(NT) IS 4’ 
GOTO 999 
ENDIF 
READ(2,*)NPHI 
IF (NPHI.GT. 200) THEN 
WRITE(6,*) MAXIMUM NUMBER OF POINTS(NPHI) IS 200° 
GOTO 999 
ENDIF 
READ (4, *)CNPHI 
IF(CNPHI.GT. 100) THEN 
WRITE(6,*) MAXIMUM NUMBER OF POINTS(CNPHI) IS 100’ 
GOTO 999 ) 
ENDIF 
C LOAD THE WEIGHTS AND ABSCISSAS IN THE VECTORS. 
DO 1 M=1,NT 
READ(1,*,END=1)XT(M) , AT(M) 
1 CONTINUE 
DO 2 M=1,NPHI 
READ(2,*,END=2)X(M) ,A(M) 
2 CONTINUE 
DO 4 M=1,CNPHI 
READ (4, *, END=4)XP(M) , AP(M) 
4 CONTINUE 
CLOSE(1) 
CLOSE(2) 
CLOSE(4) 
WRITE(6,*) ENTER THE PLOTTING INCREMENT IN DEGREES’ 
READ(5,¥*)DT 
IF(ICALC.EQ.0) THEN 
WRITE(6,*)’ENTER HIGHEST MODE’ 
READ(S,*) MODES 
NBLOCK=2*MODES+1 
NROW=NBLOCK*(2*NP-3) 
C CHECK FOR BOUNDS VIOLATION 


IF(NROW.GT.30000) THEN 
WRITE(6,*) ’BOUNDS VIOLATION FOR ARRAY CT(.)’ 
GOTO 999 
ENDIF 
ENDIF 
WRITE(6,*)’ENTER PHI (observation) IN DEGREES’ 
READ(S,*)P 
PHI=P*RAD 
IFCICALC.EQ.0) THEN 
WRITE(6,*)’ENTER THE SCAN ANGLES IN DEGREES: THETA, PHI’ 
READ(S,*) THD,PHD 
THS=THD*RAD 
PHS=PHD *RAD 
WRITE(6,*)’ENTER COMPLEX IMPEDANCE, OHMS: (Real,Imag)’ 
READ(S,*) IMP 
PESE 
Cx*x***e*e*ek*x*xkREAD CURRENT COEFFICIENTS IF ICALCHO % «ex RKKKKKEKKEKKKEE KX 
WRITE(6,*) ’ENTER LETTER CODE FOR FILE curcoefsdat’ 
READ(S,*) ITH 
CC=’curcoefsdat’//ITH 
OPEN(31,FILE=CC) 
WRITE(6,*) ’FILE OPENED IS ’,CC 
READ(31,*) ITH,SELECTION,BASE,NP,MT,MP,N 
READ(31,*) NROW,MODES,NBLOCK,THS,PHS,ARAD,PHIO 
READC31,*) IMP,ZP,RS 
READ(31,*) (RH(L),L=1,NP) 
READ(31,*) (ZH(L),L=1, NP) 
READ(31,*) (ZLO(L),L=1,NP) 
READ(31,*) (CT(L),L=1,NROW) 
GCE@SE(31) 
WRITE(6,*) ’RUN CODE READ FROM DISC : ’,ITH 
ENDIF 
C*¥**¥***k*k*k*xkREADY TO CALCULATE THE PATTERN* x **&* * XH KKHKKKEKKKK KEKE KES 
IT=INT((STOP-START)/DT)+1 
MHI=MODES+1 
OPEN(8, FILE=OUT) 
WRITE(8,8000) 2.*BASE,NP,MODES,RS,ZP 
8000 FORMAT(//,’*** BOR RADIATION PATTERN FOR A CIRCULAR’ 
*’ DISC RADIATOR USING GENEX ***’, 
*//,2X,’BOR DIAMETER (WAVELENGTHS)=’ ,F5.2,/,2X, 
* ’NUMBER OF GENERATING POINTS (NP)=’ ,14,/,2X, 
* ’ HIGHEST MODE NUMBER USED (MODES)=’,13, 
See SUR PAGCE RADIUS ,FS.2,/7 , 2h cPRiMe. 4Fo. 2) 
WRITE(8,30) NT,NPHI,CNPHI 
30 FORMAT(/,12X,’ NT NPHI CNPHI’, 
* PON, 15,2), 15,2%,10) 
TPC ISEG.EQ.0) WRITECS, 1300) 
1300 FORMAT(/,10X,’ INDEX’ ,8X,’Z(1I)’,10X, ’’RHO(I)’ ,12X,’SURF IMPED’ ) 
MP=NP-1 
MT=MP~-1 
N=MT+MP 
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DO 52 I=1,NP 

TFCABS(ZE(1)) LT @eot Zhe cr 

IF(ABS(RH(I)).LT..001) RH(I)=0. 

ZHB=ZH (I) /BK 

RHB=RH(1)/BK 
C ASSIGN SURFACE IMPEDANCE AT THIS POINT. THE SURFACE IMPEDANCE OF SEGMENT 
CG I 1S Zkom) 

IFCICALC.EQ.0) Z1O0(1I)=IMP/(120.*PI) 

IFCISEG.EQ.0) WRITE(8,8004) I,ZHB,RHB,IMP 
eZ CONTINUE 
8004 FORMAT(11X,14,4X,F8.3,8X,F8.3,6X,2F8.2) 
COIR ORR IOI OOOO OR IORI IOI III IG IOIOIIOIGII IO IOI IOIOIOR AOR IGG i ka ate ak 
C BIG LOOP * > * * * * * * * * * * * € OK KOK KK K KOK KK KOK OK OK OK 
C 
C MODE LOOP TO CALCULATE THE CURRENT COEFFICIENTS. POS AND NEG MODES * 

SCONE IN THE SAME ITERATION OF THE LOOP 

Cr Vv 
C 400 CONTINUE < * * * * * *© & * © *€ *& © KK € KH K€ OK KK KOK KOKO OK KR OK OX 
Ce It ok ok ee ee KZT LOAD , ZMAT , GENEX , DECOMP , SOLVE 

IFCICALC.EQ.0) THEN 

TF (CABS (CIMP) -NE.O) CALL ZLOAD(CNP. RH. ZE-ZLO. Zk 

DO 400 M=1,MHI 

NM=M-1 

CALL ZMAT(NM,NM,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 
Ce KE KEE KEKE KKK KEKE KEKE KEKE KEREKEKEEEEEEPOSTTIVE MODE +n 


e or ce t,phi_ | 
C ) @Z Z | 
C | +n +n | 
meee S| | 
G [Pe phist phi,phi | 
E [eZ Z | 
© | +n +n | 


CLES SELES E SESE SES SESS ES LSS SS SLES LES ESET ESE TESS TESS EST SSE LES TT ETT SST TS SS SF 


C MODE SYMMETRY IN n FOR ZMAT EXISTS BUI IS NOT USED INSTEES VERS Lane 
C MODE SYMMETRY IN I USED DIRECTLY FOR PREIS UF 059051605 inb (ce 


CoRR OOO OR ORR OR ROR IOI IO IOIGIG IG OR IOIOIGIOIGO a HOR R (Ok a tGK 


CG - | wigeee eee le = et t,phi | 
Cc )6hU6U C2 Z ale =o | 
G-  lteen -n | | +n +n | 
oar | = | | 
C> | ephaet phi,phi | eee phat phi,phi | 
C ieee? Z | | <2 Z 
cy ea =n | | +n +n | 
Ck ROR ROR ORR RR ERK KE KK KE EEKEKCURRENT COEFFICIENTS 
GC — SON 0 RE a ee =1 

C ee heel | | | | ean 
C [-i .|_ SiR S2> Sele > OR On| |v | 
C | +n ieeeroe| | | | +n | 
C | | (Reena poeta (252226 | | 
C ee | | | | Amy st 
C CT =| to =|) Ore Zl” OM 


0 


——_ oe ee oe —_—— = ee oe oe 


| | | 
| | | 
| | Bead! 
| vs el 
=m | | | | 
| | i a Mss ie ae | | | 
CoRR IOI OR OO ect POSITIVE MODE +n 
TR(CABSCIME) NE, O) CALL ZIOTCMT MP, ZE,Z) 
CALL GENEX(NM,NP,NT,NPHI,CNPHI,XT,AT,X,A, 
* XP,AP,THS,PHS,ARAD,RH,ZH,B) 
CALL DECOMP(N,IPS,Z) 
CAEL SOLVE(N. IPS ;Z,B,¢) 
NTOP1=MODES-NM 
NTOP2=NBLOCK-(NTOP1+1) 
NS2=NTOP1*N 
NS 1=NTOP2*N 
NROW=NBLOCK*N 
OSES SESE SESS ST SSS et ee eee ee eS SSE SSCS SSCS SP Pe SES TST SST SS 
C STORE CURRENT COEFFICIENTS IN ONE LONG VECTOR CT(.). MOST NEGATIVE 
CeMODE IS AT THE TOP; MOST POSITIVE AT THE BOTTOM. UNIT 30 IS FORMATED 
C OUTPUT; UNIT 31 IS FREE FORMATED (USED BY GAIN PROGRAM). 
Ck RRR OOO RR OO a iii oii toi ic ak 
OPEN (30 ,FILE=CURCO) 
WRITE(30,5027)’Mode(’,NM,’) = ’,NM,’NM = ’,NM 
DO 401 L=1,N 
WRITE OR S026) °CTOSLt NS1,7) = ,C(L) 
401 CT(L+NS1)=C(L) 
5027 FORMAT(/A,1I,A,I/A,I/) 
5026 FORMAT(5X,A,14,A,F10.4,’+ (’,F10.4,’)*i’) 
Cee KEKE KKX*EEKEND OF POSITIVE MODE*******BEGIN NEGATIVE MODE -n 
C MODE SYMMETRY IS USED IF ISYM=0 
IF(NM.NE.O) THEN 
NMN=-NM 
C USE BRUTE FORCE IF ISYM # 0 
IFCISYM.NE.0O) THEN 
CALL GENEX(NMN,NP,NT,NPHI,CNPHI,XT,AT,X,A, 
* XP,AP,THS,PHS,ARAD, RH, ZH,B) 
CALL ZMATCNMN,NMN,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 
PeaCGABSCIMP) = NE.O) CALL ZTOT(MT,MP,ZL,Z) 
CALL DECOMP(N,IPS,Z) 
Cea, SOLVECN, 1PS;,2,B;C) 
BENDIS 
WRITE( 30,5027) ’Mode(’,NMN,’) = ’,NMN,’NM = °,NM 
DO 402 L=1,MT 
CT(L+NS2)=+C(L) 


=T 


AAAAANA 


WRITE(30,5026) ’CT(’,L+NS2,’) = ’,CT(L+NS2) 
402 CONTINUE 

DO 403 L=MT+1,N 

IF(ISYM.NE.0O) CT(L+NS2)=+C(L) 

IFCISYM.EQ.0) CT(L+NS2)=-C(L) 

WRITE(30,5026) ’CT(’,L+NS2,’) = ’,CT(L+NS2) 


i. 


403 CONTINUE 
CER KK KEE EEK REE KEE KKK ERK EK EKER KEKEKEKERKEKEEND OF NEGATIVE MODE -nr 
ENDIF 
400 CONTINUE 
CLOSE(30) 
Co III RIO RIOR IO RII III IO ORR III IG OR COI ICI 
C WRITE THE VECTOR OF CURRENT COEFFICIENTS ON DISC FOR GAIN CALC. 
C (FREE FORMAT TO UNIT 31.) 
CHEEK KKK KKK KE EKE KKK KEE KKK KEK EKKEEEKK EK KEKE KE EKE KKK KEK KEEK EE KEK EEK EK EKEKKKKEKEK 
OPEN(31,FILE=CC) 
WRITE(31,*) ITH,SELECTION,BASE,NP,MT,MP,N 
WRITE(31,*) NROW,MODES ,NBLOCK,THS,PHS,ARAD, PHI 
WRITE(31,*) IMP,ZP,RS,IFF 
WRITE(31,*) (RH(L),L=1,NP) 
WRITE(31,*) (ZH(L) ,L=1,NP) 
WRITE(31,*) (ZLO(L) ,L=1,NP) 
WRITE(31,*) (CT(L),L=1,NROW) 
WRITE(31,*) GNT,GNPHI,CGP 
WRITE(6,*) *’CURRENT COEFFICIENTS WRITTEN ON DISC’ 
CLOSE(31) 
ENDIF 
CoRR IIR GIGI III IG III III III IR IO III IG a CICK 
C END OF BIG LOOP x* > * * eK Ke KK KK KK KK KK KE K K KK OK K OK OK OK OK 
ie Vv 
Cx * * kK Ke eC * * * * kK K€ KK KOR OF * * * kK K eK KK KK K KK K 
CREEK KKK EK KKK KK EKKK KEK KEK KEKE KEKE EERE EEK KEKE KK KEKKKKERK KKK KKKKKEKEKEKEK 
G 
C BEGIN PATTERN LOOP FROM START TO STOP IN INCREMENTS OF DT (ALL IN DEG) 
G 
DO 500 [=iyr 
THETA=START + DT*FLOAT(I-1) 
THX=THETA*RAD 
PHR=PHI 
IF(THETA.GT.180.) THEN 
PHR=PHI+PI 
THX=(360.-THETA) *RAD 
ENDIF 
Er 1=(0n, 0 
EP1=(0.,0 
ETZ= (Om or 
EP2=(0. ,0 
DO 300 M= 
NM=M-1 
EXP1=CEXP(CMPLX(0O. , FLOAT(NM) *PHR) ) 
EXP2=CONJG(EXP1) 
CoRR ORR OKO OK ROK fOr ror rok ko kok ok ok kkk ok oko PL ANE 
CALL PLANE(NM,NM,NP,NT,RH,ZH,XT,AT,THX,R) 
NTOP1=MODES-NM 
NTOP2=NBLOCK-(NTOP1+1) 
NS2=NTOP1*N 
NS1=NTOP2*N 


—~] 
tO 


DOm250 L=1. MT 
ET1=ET1+R(L) *CT(L+NS1) *EXP1i 
EP1=EP1+R(L+N) *CT(L+NS1) *EXP1 
IF(NM.EQ.0) GOTO 250 
ET1=ET1+R(L) *CT(L+NS2) *EXP2 
EP1=EP1-R(L+N) *CT(L+NS2) *EXP2 
250 CONTINUE 
DO 260 L=1,MP 
ET2=ET2+R(L+MT) *CT(L+NS1+MT) *EXP1 
EP2=EP2+R(L+MT+N ) *CT(L+NS1+MT) *EXP1 
IF(NM.EQ.0) GO TO 260 
ET2=ET2-R(L+MT) *CT(L+NS2+MT) *EXP2 
EP2=EP2+R(L+MT+N) *CT(L+NS2+MT ) *EXP2 
260 CONTINUE 
300 CONTINUE 
C FEED CONTRIBUTION IN THE FAR FIELD IS ETF,EPF 
CREEK KHER KERR ERE EKER KEKE KEKE KEK AK AKKEKEKERESS J 1 , 
C USE CLOSED FORM FEED EXPRESSION FROM SUBROUTINE ANTFF IF 
C IFF=0; OTHERWISE USE BRUTE FORCE EVALUATION FROM CIRCRTP 
€ 
ROBS=1000. 
IF(IFF.EQ.0) THEN 
CALL ANTFF(THX,PHR,THS,PHS,ARAD,ETF,EPF) 
ELSE 
RHB=ROBS*SIN(THX) 
ZHB=ROBS*COS (THX) 
CALL CIRCRTP(CNPHI,XP,AP,ARAD,THS,PHS, 
* PHR,RHB,ZHB,CIRCR,CIRCT,CIRCP) 
C REMOVE THE 1/R DEPENDENCE BECAUSE EXP(-jkR)/R IS OMITTED IN 
C THE SCATTERED FIELDS ET AND EP 
ETF=CIRCT*ROBS 
EPF=CIRCP*ROBS 
ENDIF 
CL PTL SLES E SS SSS ESTES SSE SS SSS SSE STL ESS STE TS ST SS STS SS SS SF 
ETHF (I) =CABS(ETF) 
EPHF (I) =CABS(EPF ) 
ET=(ET1+ET2) *UC 
EP=(EP1+EP2) *UC 
ETSCAT(1I)=CABS (ET) 
EPSCAT(1I)=CABS(EP) 
C TOTAL E-THETA AND E-PHI COMPONENTS 
EC=ET+ETF 
EX=EP+EPF 
ECV(I)=CABS(EC) 
EXV(1I)=CABS (EX) 
ECR=REAL(EC) 
ECI=AIMAG(EC) 
EXR=REAL(EX) 
EXI=AIMAG(EX) 
PHC(I)=ATAN2(ECI ,ECR+1.e-10)/RAD 
PHX(I)=ATAN2(EXI ,EXR+1.e-10)/RAD 


id 


ANG(1)=THETA 
ECX=AMAX1(ECX, ECV(I) , EXV(I)) 
500 CONTINUE 
WRITE(6,*) ’MAX E VALUE=’ ,ECX 
WRITE(8,103) P,THS/RAD,PHS/RAD, ARAD,ECX 
103 FORMAT(/,10X,’PHI OF RECEIVER (DEG)=’,F8.2,/,10X, 


* "ANTENNA SCAN: THETA (DEG)=’,F8.2,/,10X, 
* PHI CDEG) = 25 EiSe2 a / Oke 
* "ANTENNA RADIUS: ARAD=" FS. oe7 some 
* ’MAXIMUM FIELD VALUE (V/M)=’,E15.5) 


IF(IFF.EQ.0) THEN 
WRITE(8,113) 
ELSE 
WRITE(8,213) ROBS 
ENDIF 
113 FORMAT(/,10X,’CLOSED FORM FAR FIELD PATTERN FROM ANTFF IS USED’) 
213 FORMAT(/,10X,’FAR FIELD COMPUTED USING CIRCRTP, ROBS=’ ,F10.2) 
C STORE DATA FOR NORMALIZED CO- AND CROS-POLARIZED PATTERNS 
DO 600 I=1,IT 
ECP(I)=(ECV(I) /ECX) **2 
EXP(I)=(EXV(1)/ECX) **2 
ECP(I)=AMAX1(ECP(I) ,1.E-6) 
EXP(I)=AMAX1(EXP(I),1.E-6) 
ECP(I)=10. *ALOG10(ECP(I)) 
EXP(I)=10. *ALOG10(EXP(I) ) 
600 CONTINUE 
IF(IPRINT.EQ.O) THEN 
WRITE(8,5015) 
5015 FORMAT(//,7X,’ ANGLE’ ,15X, ’CO-POL’ ,25X, °X-POL’ ,/,7X, 
1’(DEG)’ ,4xX,2(’ (VOLTS) , 4x, (DEG)”’ . 3X, "(DB=REL) 54) 
DO 9000 L=1,IT 
WRITE(8,5016) ANG(L),ECV(L) ,PHC(L),ECP(L) ,EXV(L) , PHX(L) 
1,EXP(L) 
5016 FORMAT(5X,F6.2,3X,2(F8.4,3X,F7.2,3X,F7.2,3X)) 
9000 CONTINUE 
ENDIF 
OPEN(2,file=AN) 
OPEN(3,file=CPOL) 
OPEN(4,file=XPOL) 
OPEN(7,file=ETSCATTER) 
OPEN(8 ,file=EPSCATTER) 
OPEN(9,file=ETFEED) 
OPEN(10,file=EPFEED ) 
DO 9097 I=1,1T 
WRITE(2,5019)I, ANG(I) 
WRITE(3,5020)I, ECP(I) 
WRITE(4,5021)I, EXP(TI) 
WRITE(7 ,5022)1I, ETSCAT(I) 
WRITE(8,5023)I, EPSCAT(I) 
WRITE(9,5024)I, ETHF(I) 
9097 WRITE(10,5025)I, EPHF(I) 


PGi FORMATG ANG(?,1,*)=')F8.3,°;") 
SO20 PORMAI ECP(’ ,1;’)=" F8.3,7;?) 
Bedi FURMATY EXP(’,1,’)=",F8.3,°;*) 
BO22 FORMAT ( EISCAT(’ ,1,’)=',F8.3,°:*) 
5023 FORMAT(’EPSCAT(’ ,I,’)=',F8.3,°:") 
5024 FORMAT(’ETHF(’,I,’)=’,F8.3,’;’) 
Some MORMAn EPH (*.1,?)=?,F8.3,.°;?) 
CEGSEC2) 
CROSE C3) 
CLOSE(4) 
CHOSE (7) 
CLOSE(8) 
CLOSE(9) 
CLOSE(10) 
999 CONTINUE 
STOP 
END 
Ce RRR RRR RR KKK KKK RR RE KK KK RK KEK KK EK KK KERR KK KKEKEKKKEKKKEND OF MAIN PROGRAM. 
C 
Ce ee KKK RK KK ROKK RK KKK KEK KK KKK KEKEKEKEKKSUBROUTINE ZMAT. 


C REFERENCE: AN IMPROVED E-FIELD SOLUTION FOR CONDUCTING BOR 


C J.R.MAUTZ AND R.F. HARRINGTON 
C TECHNICAL REPORT TR-80-1 

C ROME AIR DEVELOPMENT CENTER 

C CONTRACT NO. F-30602-79-C-0011 


SUBROUTINE ZMAT(M1,M2,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 


COMPUTE THE MM IMPEDANCE MATRIX ELEMENTS. THIS IS FROM MAUTZ AND 
HARRINGTON (NO CHANGES). 


sO) 0} °C) 


COMPLEX Z(100000) ,U1,U2,U3,U4,U5,U6,U7,U8,U9,UA,UB,G4A(4) ,G5A(4) 
COMPLEX G6A(4) ,G4B(4) ,GS5B(4),G6B(4) ,H4A,HS5A,H6A,H4B,H5B 
COMPLEX H6B,UC,UD,GA(100) ,GB(100) 
DIMENSION RH(400) ,ZH(400) ,X(100) ,A( 100) ,XT(4) , AT(4) 
DIMENSION RS(400) ,ZS(400) ,D(400) ,DR( 400) ,DZ(400) 
DIMENSION DM(400) ,C2(100) ,C3(100) ,R2(4) ,Z2(4) 
DIMENSION C4(100) ,C5(100) ,C6(100) ,27(4) ,R7(4) ,Z8(4) ,R8(4) 
ero. 
CP=.1 
DO 10 I=2,NP 
I2=I-1 
RS(I2)=.5*(RH(I)+RH(I2) ) 
ZS(I2)=.5*(ZH(1)+ZH (12) ) 
Di=.5*(RH(I)-RH(I2)) 
D2=.5*(ZH(I)-ZH(I2)) 
D(I2)=SQRT(D1*D1+D2*D2) 
DR(I2)=D1 
DZ(I2)=D2 
Wr chs G12). 600.) RS(12)=1. 
Dail?) =DtT2) /RS (12) 
10 CONTINUE 


Gt 


29 
Ti 


M3=M2-M1+1 
M4=M1-1 
PI2=1.570796 

DO 11 K=1,NPHI 
PH=PI2*(X(K) +1.) 
C2(K)=PH*PH 
SN=SIN(.5*PH) 
C3(K)=4.*SN*SN 
A1=PI2*A(K) 
D4=.5*A1*C3(K) 
DS5=A1*COS(PH) 
D6=A1*SIN(PH) 
M5=K 

DO 29 M=1,M3 
PHM=(M4+M) *PH 
A2=COS (PHM) 
C4(M5S)=D4*A2 
C5(MS)=D5*A2 


C6(M5 ) =D6*SIN(PHM) 


MS=MS+NPHI 
CONTINUE 
CONTINUE 

MP=NP-1 

MT=MP-1 

N=MT+MP 

N2N=MT*N 

N2=N*N 

Ul=(Orm aso) 

U2=(0 2.) 
JN=-1-N 

DO 15 JQ=1,MP 
KO-2 

IF(JQ.EQ.1) KQ=1 
IF(JQ.EQ.MP) KQ=3 
R1=RS(JQ) 
Z1=ZS(JQ) 
D1=D(JQ) 

D2=DR( JQ) 

D3=DZ( JQ) 
D4=D2/R1 

DS=DM( JQ) 
SV=D2/D1 
CV=D3/D1 
T6=CT*D1 
T62=T6+D1 
T62=T62*T62 
R6=CP*R1 
R62=R6*R6 

DO 2) =i NT 
R2(L) =R1+D2*XT(L) 
Z2(L)=Z1+D3*XT(L) 


i? 


26 


at 
28 
42 


40 


3D 


34 


CONTINUE 
U3=D2*U1 

U4=D3*U1 

DO 16 IP=1,MP 
R3=RS(IP) 

Z3=ZS(IP) 

R4=R1-R3 

Z4=Z1-Z3 
FM=R4*SV+Z4*CV 

PHM=ABS (FM) 
PH=ABS(R4*CV-Z4*SV) 
D6=PH 

IF(PHM.LE.D1) GO TO 26 
D6=PHM-D1 
D6=SQRT(D6*D6+PH*PH ) 
IF(IP.EQ.JQ) GO TO 27 
Kp=t 

IF(T6.GT.D6) KP=2 
IF(R6.GT.D6) KP=3 


GU tO) 28 

KP=4 

GO TO (41,42,41,42) ,KP 
DO 40 L=1,NT 
D7=R2(L)-R3 
D8=Z2(L)-Z3 


Z7 (L)=D7*D7+D8*D8 
R7(L)=R3*R2(L) 

ZS (ia) =225*77(L) 
R8(L)=.25*R7(L) 
CONTINUE 
Z4=R4*R4+Z4*24 
R4=R3*R1 

RS=.5*R3*SV 

DO 33 K=1,NPHI 
A1=C3(K) 

RR=Z4+R4*A1 

UA=0. 

UB=0. 

IF(RR.LT.T62) GO TO 34 
DO 35 L=1,NT 
R=SQRT(Z7(L)+R7(L)*A1) 
SN=-SIN(R) 

CS=COS(R) 
UC=AT(L)/R*CMPLX(CS,SN) 
UA=UA+UC 
UB=XT(L)*UC+UB 
CONTINUE 

GO TO 36 

Does? e-1, NT 
R=SQRT(Z8(L)+R8(L) *A1) 
SN=-SIN(R) 


=] 
—~J 


Si 


38 
BS, 


36 


33 


46 


45 


CS=COS(R) 
UC=AT(L)/R*SN*CMPLX(-SN,CS) 
UA=UA+UC 
UB=XT(L)*UC+UB 
CONTINUE 

A2=FM+RS*A1 
D9=RR-A2*A2 
R=ABS(A2) 

D7=R-D1 

D8=R+D1 

D6=SQRT (D8*D8+D9) 
R=SQRT(D7*D7+D9) 
IF(D7.GE.0.) GO TO 38 
A1=ALOG((D8+D6) *(-D7+R)/D9)/D1 
GO TO 39 
A1=ALOG((D8+D6)/(D7+R))/D1 
UA=A1+UA 
UB=A2*(4./(D6+R)-A1)/D1+UB 
GACK)=UA 

GB(K)=UB 

CONTINUE 

K1=0 

DO 45 M=1,M3 

H4A=0. 

HSA=0. 

H6A=0. 

H4B=0. 

HSB=0. 

H6B=0. 

DO 46 K=1,NPHI 
Ki=Kit+i 

D6=C4(K1) 

D7=CS(K1) 

D8=C6(K1) 

UA=GA(K) 

UB=GB(K) 
H4A=D6*UA+H4A 
HSA=D7*UA+HSA 
H6A=D8*UA+H6A 
H4B=D6*UB+H4B 
HSB=D7*UB+HSB 
H6B=D8*UB+H6B 
CONTINUE 

G4A(M)=H4A 
GSA(M)=HSA 
G6A(M)=H6A 
G4B(M)=H4B 
GSB(M)=HSB 
G6B(M)=H6B 

CONTINUE 

IF(KP.NE.4) GO TO 47 


=] 
CO 


65 


64 


66 
63 


67 


41 


25 


Ls 


62 


§1 


A2=D1/(PI2*R1) 

D6=2./D1 

D8=0. 

DO 63 K=1,NPHI 
A1=R4*C2(K) 

R=R4*C3(K) 

IF(R.LT.T62) GO TO 64 
D0. 

DO 65 L=1,NT 
D7=D7+AT(L)/SQRT(Z7(L)+A1) 
CONTINUE 

GO TO 66 
A1=A2/(X(K) +1.) 
D7=D6*ALOG(A1+SQRT(1.+A1*A1)) 
D8=D8+A(K)*D7 

CONTINUE 

A1=.5*A2 

A2=1./A1 
D8=-PI2*D8+2./R1*(BLOG(A2) +A2*BLOG(A1) ) 
DO 67 M=1,M3 
GSA(M)=D8+G5A(M) 
CONTINUE 

GO TO 47 

DO 25 M=1,M3 

G4A(M)=0. 

GSA(M)=0. 

G6A(M)=0. 

G4B(M)=0. 

GSB(M)=0. 

G6B(M)=0. 

CONTINUE 

DOmtoeL=1.NT 

A1=R2(L) 

R4=A1-R3 

Z4=Z2(L)-Z3 
Z4=R4*R4+Z4*24 

R4=R3*A1 

DO 17 K=1,NPHI 
R=SQRT(Z4+R4*C3(K)) 
SN=-SIN(R) 

CS=COS(R) 
GA(K)=CMPLX(CS,SN)/R 
CONTINUE 

D6=0. 

IF(R62.LE.Z4) GO TO 51 
DO 62 K=1,NPHI 
D6=D6+A(K)/SQRT(Z4+R4*C2(K) ) 
CONTINUE 
Z4=3.141593/SQRT(Z4/R4) 
D6=-PI2*D6+ALOG(Z4+SQRT(1.+Z4*Z4) )/SOQRT(R4) 
A1=AT(L) 


SZ 


30 
is 
47 


A2=XT(L)*A1 

K1=0 

DO 30 M=1,M3 

US=0:. 

U6=0. 

Uz=0. 

DO 32 K=1,NPHI 
UA=GA(K) 

KiKi 1 

U5=C4(K1) *UA+U5 
U6=C5 (K1) *UA+U6 
U7=C6(K1) *UA+U7 
CONTINUE 

U6=D6+U6 
G4A(M)=A1*U5+G4A(M) 
GSA(M) =A1*U6+G5A(M) 
G6A(M)=A1*U7+G6A(M) 
G4B(M) =A2*U5+G4B(M) 
G5B(M)=A2*U6+G5B(M) 
G6B(M)=A2*U7+G6B(M) 
CONTINUE 

CONTINUE 

A1=DR(IP) 

UA=A1*U3 
UB=DZ( IP) *U4 
A2=D(IP) 

D6=-A2*D2 

D7=D1*A1 . 
D8=D1*A2 

JM=IN 

DO 31 M=1,M3 
FM=M4+M 
A1=FM*DM(IP) 
H5A=GS5SA(M) 
HSB=G5B(M) 
H4A=G4A(M)+H5A 
H4B=G4B(M)+H5B 
H6A=G6A(M) 
H6B=G6B(M) 
U7=UA*HSA+UB*H4A 
US8=UA*HSB+UB*H4B 
U5=U7-U8 

U6=U7+U8 
U7=-U1*H4A 
U8=D6*H6A 
U9=D6*H6B-A1¥*H4A 
UC=D7 * (H6A+D4*H6B) 
UD=FM*D5*H4A 
K1=IP+JM 

K2=Ki+1 

K3=Ki+N 


K4=K2+N 
KS=K2+MT 
K6=K4+MT 
K7=K3+N2N 
K8=K4+N2N 
GO TO (18,20,19),KQ 
18 Z(K6)=U8+U9 
TCP bO.1) Ga TO, 21 
Z(K3)=Z(K3)+U6-U7 
Z(K7 )=Z(K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 
21 Z(K4)=U6+U7 
Z(K8) =UC+UD 
GO TO 22 
19 Z(K5)=Z(KS)+U8-U9 
IF(IP.EQ.1) GO TO 23 
Z(Ki)=Z (Ki )+US+U7 
Z(K7 )=Z(K7 )+UC-UD 
IF(IP.EQ.MP) GO TO 22 
23 Z(K2)=Z(K2)+U5-U7 
Z(K8 )=UC+UD 
Go 10) 22 
20 Z(K5)=Z(K5)+U8-U9 
Z(K6 ) =U8+U9 
IF(IP.EQ.1) GO TO 24 
Z(K1)=Z(K1i)+US+U7 
Z(K3)=Z(K3) +U6-U7 
Z(K7 )=Z (K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 
24 Z(K2)=Z(K2)+US-U7 
Z(K4) =U6+U7 
Z(K8)=UC+UD 
22 Z(K8+MT )=U2*(D8* (H5A+D4*H5B)-A1*UD) 
JM=JM+N2 
31 CONTINUE 
16 CONTINUE 
IN=IN+N 
15 CONTINUE 
RETURN 
END 
Ce KEKE EK KEKE EEK EEE KEKE KE KEKE EKEESEKESEESUIRROUTINE SOLVE. 
C REFERENCE : TECHNICAL REPORT TR-80-1 
SUBROUTINE SOLVE(N,IPS,UL,B,X) 


C SEE MAUTZ & HARRINGTON FOR DETAILS 


COMPLEX UL(100000) ,B(800) ,X(800) ,SUM 
DIMENSION IPS(800) 

NP=N+1 

IP=IPS(1) 

X(1)=BCIP) 


DO 2 I=2,N 
IP=1PsS( 1) 
IPB=IP 
IMi=I-1 
SUM=(0.,0.) 
DOwe J =1,.061 
SUM=SUM+UL (IP) #X(J) 
1 IP=IP+N 
2 X(I)=BCIPB)-SUM 
K2=N*(N-1) 
IP=IPS(N)+K2 
X(N)=X(N)/ULCIP) 
DO 4 IBACK=2,N 
I=NP-IBACK 
K2=K2-N 
IPI=IPS(1I)+K2 
PPi- iri 
SUM=(0. ,0.) 
IP=IPI 
DO 3 J=IPi1,N 
IP=IP+N 
3 SUM=SUM+UL( IP) *X(J) 
4 X(I)=(X(1)-SUM)/ULCIPI) 
RETURN 
END 
CREEK KEKE HEEEKEESRESEEEEEXSIIBROUTINE DECOMP 
C REFERENCE: TECHNICAL REPORT TR-80-1 
SUBROUTINE DECOMP(N,IPS,UL) 
C 
C SEE MAUTZ & HARRINGTON FOR DETAILS 
c 
COMPLEX UL(100000) , PIVOT, EM 
DIMENSION SCL(800) , IPS(800) 
DO & I=1,N 
IPS(I)=I 
RN=0. 
Ji=I 
DO 2 J=1,N 
ULM=ABS (REAL(UL(J1)))+ABS(CAIMAG(UL(J1))) 
Ji=Ji+N 
IF(RN-ULM) 1,2,2 
1 RN=ULM 
2 CONTINUE 
SCL(I)=1./RN 
5 CONTINUE 
NMi=N-1i1 
K2=0 
DO 17 K=1,NM1 
BIG=0. 
DO 11 I=K,N 
IP=IPS(I) 


GO 
to 


10 


Gl 


14 


15 


18 


16 


Vi 


IPK=IP+K2 
SIZE=(ABS(REAL(UL(IPK) ) )+ABS(AIMAG(UL(IPK))))*SCL(IP) 
IF(SIZE-BIG) 11,11,10 
BIG=SIZE 

IPV=I 

CONTINUE 

IF(IPV-K) 14,15,14 
J=IPS(K) 
IPS(K)=IPS(IPV) 
IPS(IPV)=J 
KPP=IPS(K)+K2 
PIVOT=UL(KPP) 

KP1=K+1 

DO 16 I=KP1,N 

KP=KPP 

IP=IPS(I)+K2 
EM=-UL(IP)/PIVOT 
UL(IP)=-EM 

DO 16 J=KP1,N 

IP=IP+N 

KP=KP+N 

UL( IP) =UL(IP)+EM*UL (KP) 
CONTINUE 

K2=K2+N 

CONTINUE 

RETURN 

END 


Cee KERR KEKEKKEPUNCTION BLOG 


C REFERENCE: TECHNICAL REPORT TR-80-1 ;(page 56) 


FUNCTION BLOG(X) 

Pe Cx.Gl 1) co. To 1 

X2=X*X 

BLOG=( ( .075*X2-. 1666667) *X2+1. ) *X 
RETURN 

BLOG=ALOG(X+SQRT(1.+X*X) ) 

RETURN 

END 


CHEER ERE EEE EERE EEE EERE RARER EE EREREREE KEKE KEKSKEKSUBROUTINE PLANE 
C REFERENCE: TECHNICAL REPORT TR-80-1 ; (pages 57-64) 


AaAaaAaANna 


SUBROUTINE PLANE(M1,M2,NP,NT,RH,ZH,XT,AT,THR,R) 


PLANE WAVE EXCITATION VECTOR IN THE FAR FIELD. FROM MAUTZ AND HARRINGTON. 
MODIFIED TO DO ONLY ONE ANGLE AND FREQUENCY PER CALL. 
x* NEW BESSEL FUNCTION ROUTINE FROM NUM. RECIPES ** 


COMPLEX R(1600),U,U1,UA,UB,FA(1500) ,FB(1500) ,F2A,F2B,F1A,F1B 
COMPLEX U2,U3,U4,U5 

DIMENSION RH(400) ,ZH(400) ,XT(4) ,AT(4) ,R2(4) ,22(4) 

MP=NP-1 

MT=MP-1 

N=MT+MP 


i3 


Poe 


N2=24N 
CC=COS(THR) 

SS=SIN (THR) 
U=(0.,1.) 
U1=3.141593*U**M1 
M3=Mi+1 

M4=M2+3 

IF(M1.EQ.0) M3=2 
MS=M1i+2 

M6=M2+2 

DO 12 IP=1,MP 

K2=IP 

I=IP+1 
DR=.5*(RH(I)-RHCIP) ) 
DZ=.5*(ZH(1)-ZH(IP) ) 
D1=SQRT(DR*DR+DZ*DZ) 
Ri=.25*(RH(1)+RH(IP)) 
IF(R1.EQ.0.) Ri=1. 
Z1=.5*(ZH(1)+*ZHCIP)) 
DR=.5*DR 

D2=DR/R1 

DOs. L=1-NT 
R2(L)=R1+DR*XT(L) 
ZO =Z1+DZ¥X TL) 
CONTINUE 

D3=DR*CC 

D4=-DZ*SS 

DS=D1¥*CC 

DO 23 M=M3,M4 
FA(M)=0. 

FB(M)=0. 

CONTINUE 

DO 15 L=1,NT 
X=SS*R2(L) *2. 
ARG=Z2(L) *CC 
UA=AT(L)*CMPLX(COS( ARG) , SIN(ARG) ) 
UB=XT(L)*UA 


C THIS LINE REPLACES HARRINGTON’S BLOCK 


25 
is 


Zo 


DO 25 M=M3,M4 
BES=BESSJ(M-2,X) 
FA(M)=BES*UA+FA(M) 
FB(M)=BES*UB+FB(M) 
CONTINUE 

CONTINUE 
IF(M1.NE.0) GO TO 26 
FA(1)=-FA(3) 
FB(1)=-FB(3) 

UA=U1 

DO 27 M=MS,M6 
M7=M-1 

M8=M+1 


s+ 


F2A=UA*(FA(M8)+FA(M7) ) 
F2B=UA*(FB(M8)+FB(M7) ) 
UB=U*UA 
F1A=UB*(FA(M8)-FA(M7) ) 
F1B=UB*(FB(M8)-FB(M7) ) 
U4=D4*UA 
U2=D3*F1A+U4*FA(M) 
U3=D3*F1B+U4*FB(M) 
U4=DR*F2A 
US=DR*F2B 
Ki=K2-1 
K4=Ki+N 
KS=K2+N 
R(K2+MT)=-D5*(F2A+D2*F2B) 
R(K5+MT)=D1*(F1A+D2*F1B) 
IF(IP.EQ.1) GO TO 21 
R(K1)=R(K1)+U2-U3 
R(K4) =R(K4) +U4-U5 
IF(IP.EQ.MP) GO TO 22 

21 R(K2)=U2+U3 
R(KS) =U4+U5 

22 K2=K2+N2 
UA=UB 

27 CONTINUE 

12 CONTINUE 
RETURN 
END 


CHEEK KEKE KEKE KEKE EEK KEKE EKER EKER ERE EEREEEEEKEKESURBROUTINE ZLOAD. 


C REFERENCE: COMPUTATION OF RADIATION AND SCATTERING FOR 


C LOADED BODIES OF REVOLUTION 

C HARRINGTON AND MAUTZ 

C REPORT AFCRL-70-0046 

C AIRFORCE CAMBRIDGE RESEARCH LABS 
C CONTRACT NO. F-19628-68-C-0180 


SUBROUTINE ZLOAD(NP,RH,ZH,ZO,Z) 


COMPUTES IMPEDANCE MATRIX ELEMENTS FOR LOADED BODIES OF REV 
ZO(I) IS THE SURF IMPEDANCE OF THE ITH SEGMENT (NP-1 SEGMENTS) 
Z(.) ARE THE IMPEDANCE MATRIX TERMS (TRIDIAGONAL FOR T-T 
SUBMATRIX; DIAGONAL FOR P-P SUBMATRIX). STORED IN COL VECTOR. 


Vaart at OP ae ee a Jo @ 


COMPLEX C1,C2,Z0(400) ,Z(2400) ,X1,X2,X3,Y1,Y2, Y3,FN(400) 
COMPLEX U1,U2,U3,XI,YI 

DIMENSION RH(400) ,ZH(400) ,RS(400) ,D(400) ,SV(400) 
PI=3.14159 

MT=NP-2 

MP=NP-14 

N=MT+MP 

DO 10 IP=2,NP 

II=IP-1 

DR=RH(IP)-RH(IT) 


CA 
art 
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DZ=ZH(IP)-ZH(II) 
D(II)=SQRT(DR*DR+DZ*DZ) 
SV(II)=DR/DCII) 
RS(II)=.5*(RH(IP)+RH(II)) 
DS=D(II)*SV(II)/2. 
Q1=RS(II)+DS 
Q2=RS(II)-DS 
FN(II)=1. 
IF((ABS(Q2).GT.1.E-6). AND. (ABS(Q1).GT.1.E-6)) 
* FN(II)=AL0G(Q1/Q2) 
10 CONTINUE 
LO=MT*3-2 
DO 20 I=1,MP 
C1=PI*Z0(I) 
IF(I.EQ.MP) GO TO 80 
KI=2 
TECE. BGs ekr=1 
IF(I.EQ.MT) KI=3 
II=I+1 
€2-pie720(01) 
A=sv 01> 
IF(ABS(A).LT.1.E-6) GO TO 41 
X1=C1*FN(I)/2./A 
X2=C1*2./A*(1.-RS(I)*FN(I)/D(I)/A) 
X3=-X2*RS(I)/D(CI)/A 
GO TO 42 
41 CONTINUE 
X1=Ci7 2278S (1) <P> 
KO =(0 Os) 
X3=C1*D(I)/6./RS(I) 
42 CONTINUE 
A=SVCIE) 
IF(ABS(A).LT.1.E-6) GO TO 45 
Y1=C2*FN(II)/2./A 
Y2=C2*2./A*(1.-RSCII)*FN(CII)/DCII)/A) 
Y3=-Y2*8S(II)/D(1i)/A 
GO TO 40 
45 CONTINUE 
Yi=6272 7RSGie <p Cir 
Y2~(OmnOe) 
Y3=C2*D(II)/6./RS(II) 
40 CONTINUE 


EFINE TRIDIAGONAL ELEMENTS FOR T-T SUBMATRIX (STORED IN COLS) 
Ul= DIAG; U2- LOWER; Us— UPPER) 


XI=X1+X2+X3 

YI=Y1-Y2+Y3 

IF(RH(I).LT.1.E-6) XI=C1/SV(I) 
IF(RHC(II).LT.1.E-6) YI=C2/SV(ITI) 
U1=XI+YI 


Y2=%1-xX3 
U3=Y1-Y3 
L=2+(1=2)%3 
PPC KI ee Oe) eG 
bi=bF1 
L2=L+2 
L3=L+3 
go to (50,60,70),ki 
50 Z(L1)=U1 
Zeke i-U2 
GO TO 80 
60e2(L1 )=U3 
Z(L2)=U1 
Z(L3)=U2 
GO TO 80 
Woe? GL lj=US 
Z(L2)=U1 
80 Z(LO+I)=2.*C1*D(I)/RS(I) 
20 CONTINUE 
RETURN 
END 
SUBROUTINE ZTOT(MT,MP,ZL,Z) 
C 
C ADDS THE SURF IMPEDANCE TERMS TO THE TRIDIAGONAL ELEMENTS OF 
C THE BOR IMPEDANCE MATRIX Z. 
€ 
COMPLEX ZL(2400) ,Z(100000) 
N=MT+MP : 
MO=MT*3-2 
DO 100 I=1,MP 
LO=MT*N+(1I-1)*N+MT 
IF(I.EQ.MP) GO TO 80 
KI=2 
me Cineg.t) Ki=1 
IF(I.EQ.MT) KI=3 
L2=(1-1)*N+I 
Esl Sjepze sh 
L3=L2+1 
M=2+3*(I-2) 
IF(KI.EQ.1) M=0 
Mi=M+1 
M2=M+2 
M3=M+3 
go to (50,60,70),ki 
50 Z(L2)=Z(L2)+ZL(M1) 
Z(L3)=Z(L3) +ZL (M2) 
GO TO 80 
60 °Z0L1)=Z(L1)+ZL (M1) 
Pile )=Z2Z (2) +Zb4 M2) 
Z(L3)=Z(L3)+ZL(M3) 
GO TO 80 


GO 
—~] 


70 ZL =Z2Z0L 1) Zea) 
712) =Z 2 ee a) 
80 Z(LO+I)=Z(LO+I)+ZL(MO+I) 
100 CONTINUE 
RETURN 


END 
COO cee SUBROUTINE OGIVE 


SUBROUTINE : OGIVE 

DATE [ea SEPTEMBER 1991 
PROGRAMMER : R.M. FRANCIS 
REVISED : 26 JANUARY 1992 
COMMENTS 


THIS SUBROUTINE WILL GENERATE DATA FOR A BODY OF REVOLUTION (BOR) IN 
THE FORM OF AN OGIVE. 
DIMENSIONS ARE NORMALIZED TO WAVELENGTH, SEE PAGE 14 FRANCIS THESIS. 
ZH = Z CUO-ORDINATE * 2*PI 
RH = RADIUS *2*PI 

SUBROUTINE OGIVE(NP,ARAD,ZH,RH,BASE,RS,ZP) 


Ore C1 Mam es Oe a 


C NP = NUMBER OF POINTS ON THE OGIVE SURFACE, MAXIMUM = 1000 

C2ZP = ZPRIME, THE POSITION ON Z WHERE THE RADIUS OF CURVATURE STARTS 
C BASE = BASE RADIUS 

C RS = RADIUS OF CURVATURE IN THE RZ,Z PLANE OF THE OGIVE 


INTEGER I,NP,NPBASE 
REAL RH(400) ,ZH(400),ZP,BASE,RS,ZCOORD,RADIUS, AL 
PI=3.1415926 
C INPUT THE VARIABLES FOR THE OGIVE,ZP,B,RS,NP 
WRITE(6,*) ’ENTER SURFACE CURVATURE (wavelengths) ’ 
READ(5,*) RS ; 
WRITE(6,*) ’ENTER ZPRIME,WHERE CURVATURE STARTS (wavelengths)’ 
READ(5,*) ZP 
WRITE(6,*) ’ENTER BASE RADIUS (wavelengths) ’ 
READ(5,*) BASE 
C PERFORM CALCULATIONS 
AL=SQRT (2. *BASE*RS-BASE**2) 
ZMAX= AL + ZP 
ANG=ASIN(AL/RS) 
L=ANINT((RS*ANG+ZP) *5) 
NP=2*L+1 


WRITE(6,*) NUMBER OF POINTS IS: ’,NP 
DZ= ZMAX/FLOAT(NP-1) 
DO 10 I=1,NP 
ZCOORD=ZMAX- FLOAT(I-1)*DZ 
IF (ZCOORD.EQ.0O.)THEN 
ZCOORD=0.0000000001 
ENDIF 
ZH(I)=2.*PI*ZCOORD 
IF (ZCOORD.LE.ZP) THEN 
RADIUS=BASE 
ELSE 
RADIUS=SQRT (RS ** 2-(ZCOORD-ZP ) **2)+(BASE-RS) 


88 


IF (RADIUS.EQ.0.)THEN 
RADIUS=0.0000000001 
ENDIF 
ENDIF 
RH(I)=2.*PI*RADIUS 
10 CONTINUE 
RETURN 
END 


CHEEK ER ER EK EEE EEE KER EEK KEKE EKEKECEKEKE EE KEKE KESIIRBROUTINE GENEX 
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SUBROUTINE : GENEX 
DATE : 14 JANUARY 1992 
REVISED peit February 1992 
PROGRAMMER <: R.M. FRANCIS 
SUBROUTINE GENEX(MODE,NP,NT,NPHI,CNPHI,XT,AT,X,A, 
* APL AP, TRS; PES, ARAD, RHO, ZHO;R) 
SOURCE ON Z AXIS AT Z=0. THIS VERSION DOES THE ANTENNA INTEGRATION 
IN RECTANGULAR COORDINATES. CNPHI POINTS USED IN X AND Y. 


EKEKKKKKKKKEKEKKKEKEKEKEKE KEKE 


* 

t t Bh * 
ev) =<W ,E> * 
at ge na * 

* 

* 

phi phi i * 
Cay) =<W ,E> * 
ae ni * 

x 


CEE KKK EEK EEE KEEKEEEKK ES 


C 
C 


THE EXCITATION VECTOR IS COMPUTED FOR THE GIVEN R,THETA AND PHI 
COMPONENTS ON SURFACE OF BOR, SPECIFIED IN SUBROUTINE CIRCRTP. 
COMPLEX R(800),CEXP,PSI,S1,$82,83,S84,S5 
COMPLEX CIRCR,CIRCT,CIRCP 
DIMENSION RH(400) ,ZH(400) ,XT(4) , AT(4) ,X( 100) ,A(100) 
DIMENSION XP(100),AP(100) ,RHO(400) , ZHO(400) 
INTEGER CNPHI 
PI=3.141592654 
BK=2.*PI 
MP=NP-1 
MT=MP-1 
N=MT+MP 


C LIMITS ON PHI INTEGRATION 


Pl=(2.*P1-0.)/2. 
P2=(2.*PI+0.)/2. 
DO) 10. 1=1 NP 

RH(I)=RHO(I)/BK 


10 ZH(1I)=ZHO(1)/BK 


DOTCO 1P=1] . AP 


C QUANTITIES FOR THE FIRST SEGMENT (POSITIVE SLOPE) 


1=IP+1 
II=IP 


DR=RH(I)-RHCII) 
DZ=ZHC 1) - Zh 
D1=SQRT(DR*DR+DZ*DZ) 
Ri=.5*(RH(I)+RH(II)) 
Z1=.5*(ZH(1I)+ZH(II)) 
SVP1=DR/D1 
CVP1=DZ/D1 
Vi=ATAN2(DR,DZ+1.E-5) 
C QUANTITIES FOR THE SECOND SEGMENT (NEGATIVE SLOPE) 
C (SKIP IF LAST SEGMENT) 
I=IP+2 
II=IP+1 
DR=RH(1I)-RHCII) 
DZ=ZH(1I)-ZHCII) 
D2=SQRT (DR*DR+DZ*DZ) 
R2=.5*(RH(1I)+RHC(II) ) 
Z2=.5*(ZH(I)+ZH(II)) 
SVP2=DR/D2 
CVP2=DZ/D2 
V2=ATAN2(DR,DZ+1.E-5) 
C BEGIN PHI INTEGRATION: R TERMS IN S1 AND S2; THETA TERM IN S3 AND S4: 
C PHI TERM IN SS. 
Si-O) 
S2=(0.,0 
So (Oe oe» 
S4=(0. ,0 
S5=(0.,0..) 
DO 20 J=1,NPHI .- 
PH=P1*X(J)+P2 
PSI=CEXP(CMPLX(0. ,-MODE*PH))*A(J) 
IF(IP.LE.MT) THEN 
C t-CURRENT INTEGRATION FOR THE POSITIVE SLOPE 
C Gauss Quadrature 
Dols Eee 
TP=Ditie ee 7 2° 
RHB=R1i+TP*SVP1 
ZHB=Z1+TP*CVP1 
TH=ATAN2(RHB,ZHB+1.E-5) 
CC=COS(V1-TH) 
SS=SIN(V1-TH) 
CALL CIRCRTP(CNPHI,XP,AP,ARAD,THS,PHS, 
* PH,RHB,ZHB,CIRCR,CIRCT,CIRCP) 
S1=S1+AT(L)*(.5+TP/D1)*CC*#CIRCR*PSI 
S2=S2+AT(L)*(.5+TP/D1)*SS*CIRCT*PSI 
13 CONTINUE 
C t-CURRENT INTEGRATION FOR THE NEGATIVE SLOPE 
C Gauss Quadrature 

DO) 14) P=iN 

TP=D2*XT(L)/2. 

RHB=R2+TP*SVP2 

ZHB=Z2+TP*CVP2 
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TH=ATAN2(RHB,ZHB+1.E-5) 

SS=SIN(V2-TH) 

6G=€0S(V2=TH) 

CALL CIRCRTP(CNPHI,XP,AP,ARAD,THS,PHS, 

* PH,RHB,ZHB,CIRCR,CIRCT,CIRCP) 
S3=S3+AT(L)*(.5-TP/D2)*CC*CIRCR*PSI 
S4=S4+AT(L)*(.5-TP/D2)*SS*CIRCT*PSI 

14 CONTINUE 
ENDIF 

C phi-CURRENT INTEGRATION 
DoOmis, L=1 NE 
TP=D1*XT(L)/2. 
RHB=Ri+TP*SVP1 
ZHB=Z1+TP*CVP1 

TH=ATAN2(RHB, ZHB+1.E-5) 

CALL CIRCRTP(CNPHI,XP,AP,ARAD,THS,PHS, 

« PH,RHB,ZHB,CIRCR,CIRCT,CIRCP) 
S5=S5+AT(L) *CIRCP*RHB/R1*PSI 

15 CONTINUE 

20 CONTINUE 

C COMPONENTS ARE STORED IN A COLUMN VECTOR: VT(1,MT), VP(MT+1,N) 
IF(IP.LE.MT) THEN 

R(IP)=(S1+S2)*D1/2. *P1+(S3+S4)*D2/2.*P1 
ENDIF 
RC IP+MT)=S5*P1*D1/2. 

30 CONTINUE 

RETURN 

END 
CEEKKKKEEKREKKEEKEKEKEKEKREE EEK KKEKEEKEKEEEREEKEKK KEE EXSIIRROUTINE CIRCRTP 


C SUBROUTINE : CIRCRTP 

C DATE : 1 OCTOBER 1997 

ey REVISED : 26 February 1992 

C PROGRAMMER : R.M. FRANCIS 

Ce2nND REVISION : 27 October 1992 

C PROGRAMMER : Ke A. KEOPP 

e 

C COMMENTS : THIS SUBROUTINE WILL RIGOROUSLY CALCULATE THE ELECTRIC FIELD 
eeeOnk A CIRCULAR APERTURE. THE FIELD IS CALCULATED AT THE COORDINATES 

C SPECIFIED BY RH(.) AND ZH(.) THE APERTURE IS LOCATED AT Z = 0, 

C AND SCANNED TO A DIRECTION (THS,PHS). THE SUBROUTINE OGIVE IS THE 

C SOURCE OF THE GEOMETRIC DATA REQUIRED BY CIRC TO PERFORM COMPUTATIONS. 

C ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 

S 

CREEK KKKK KKK KKK KEKEK EKER KKKKEE EKER EKER EKER KEKE K KK EE 
C +ant +ant * 
C - = * 
c al . \ x1 * 
C E = -jn/(4*pi*k) /_ me oe Creare )* TAPER (rp) * 
C x=0 y=0 * 
c * 
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SUBROUTINE CIRCRTP(CNPHI,XP,AP,ARAD,THS, PHS, 
* PHI ,RZ,Z,CIRCR,CIRCT,CIRCP) 
INTEGER CNPHI 
DIMENSION XP(100) ,AP(100) 
REAL K 
COMPLEX J,JK,G1,G2,X1,CON,CC,SUMX,SUMY, SUMZ 
COMPLEX CIRCR,CIRCT,CIRCP 
PI=3.141592654 
ETA=120. «Px 
SET FIELD COMPONENTS TO ZERO IN THE REAR HEMISPHERE 
CIRCR=(0.,0.) 
CIRCT=(0.40e)) 
CIREP=(0)..0-) 
Z OF THE OBSERVATION POINT IS >= O IN THE FORWARD HEMISPHERE 
IF(Z.GE.0.) THEN 
SUMX=(0.0,0.0) 
SUMY=(0.0,0.0) 
SUMZ=(0:2040 .0) 
K=6. 283185307 
I= eOr O10) 
JK= (0.0,6.283185307) 
CON=-J/(4.*PI) 
OMIT THE FACTOR OF ETA SINCE IT IS NOT INCLUDED IN SUBROUTINE ZMAT. 
ARAD**2 IS FOR SCALING OF GAUSSIAN INTEGRATION. READJUST SCALING 
TO AGREE WITH OLD SUBROUTINE SPHERE (ARBITRARY) 
CC= CON 
OUTER LOOP: INTEGRATE OVER X, -ARAD < X < ARAD. 
INNER LOOP: INTEGRATE OVER Y, -ARAD < Y < ARAD. 
EVALUATE ANTENNA FIELD AT RZ,Z: | 
XO=RZ*COS (PHI) 
YO=RZ*SIN(PHI) 
STS=SIN(-TES) 
S=RZ/SQRT (RZ**2+Z**2) 
C=Z/SQRT(RZ**2+Z¥**2) 
INTEGRATE IN X,Y COORDINATES 
DO 50 M=1,CNPHI 
X=ARAD*XP(M) 
DO 60 N=1,CNPEHI 
Y=ARAD*XP (N) 
RP=SQRT(X**2+Y **2) 
SKIP INTEGRATION POINTS BEYOND THE ANTENNA RADIUS 
IF(RP.LE.ARAD) THEN 
PHIPRIME=ATAN2(Y ,X+1.E-10) 
R=SQRT(((RZ-RP) **2+Z**2)+(4*RZ*RP*SIN( (PHI-PHIPRIME)/2) **2) ) 
G1=(((K*R) **2)-1-(JK*R) )/R**3 
G2=(3+(3*JK*R)-(K*R) **#2)/R*e*5 
X1=JK*(RP*COS(PHIPRIME-PHS) *STS-R) 
Y1=X0-X 
X2=Y1**2 
Y2=YO-Y 
SUMX=SUMX+ (CC*AP (M) *AP(N) *(G1+(X2*G2) ) *CEXP(X1) ) * TAPER (RP ) 


OZ 


SUMY=SUMY+(CC#¥AP(M) *AP(N)*Y1*Y2*G2*CEXP(X1))*TAPER(RP) 
SUMZ=SUMZ+(CC*AP(M) *AP(N) *Y1*Z*G2*CEXP(X1))*TAPER(RP) 
ENDIF 
60 CONTINUE 
50 CONTINUE 
C ASSUME AN ELEMENT FACTOR, ELMT 
c ELMT=SQRT(ABS(C)) 
ELMT=1. 
CIRCT=(C*COS (PHI )*SUMX+C*SIN (PHI) *SUMY-S*SUMZ) *ELMT 
CIRCR=(S*COS(PHI) *SUMX+S*SIN (PHI) *SUMY+C*SUMZ) *ELMT 
CIRCP=(-SIN( PHI) *SUMX+COS (PHI) *SUMY) *ELMT 


ENDIF 
RETURN 
END 
CHEE KEREKEKEKEREKKEEEKEK KEKE KEEEEEEKEEKEKERKEXESIBROUTINE TESTSPHERE 
C SUBROUTINE : TESTSPHERE 
C DATE : 21 FEBRUARY 1992 
C PROGRAMMER R. M. FRANCIS 
C COMMENTS : GENERATES A SPHERE, WITH PLOTTED POINTS 
C AT EQUAL INTERVALS IN THETA. 


SUBROUTINE TESTSPHERE(NP ,ZH,RH,BASE,RS,ZP) 
INTEGER NP,I 
REAL ZH(400) ,RH(400) , BASE,RS,ZP,SPRAD 
PI=3.1415926 
BK=2.¥*PI 
ZP=0.0 
WRITE(6,*)’ENTER NUMBER OF POINTS (NP):’ 
READ(5,*)NP 
WRITE(6,*)’ENTER SPHERE RADIUS’ 
READ(5,*)SPRAD 
BASE=SPRAD 
RS=SPRAD 
DO 1241 I=1,NP 
ANGLE=PI*FLOAT(I-1)/(2.*FLOAT(NP-1) ) 
ZH(1I)=BK*SPRAD*COS (ANGLE) 
RH(1)=BK*SPRAD*SIN (ANGLE) 
IF(RH(I).EQ.0.) THEN 
RH(I)=0.0000000001 
ENDIF 
1241 CONTINUE 
RETURN 
END 
CHEEK RK KEKE REE KEKE KKKEEEEEE EEE EEE EREREREESKEXSUBROUTINE CONE 
C GENERATES THE BOR GEOMETRY FOR A CONE. 
SUBROUTINE CONE(NP,ZH,RH,BASE,RS,ZP) 
REAL ZH(400) ,RH(400) ,BASE,RS,ZP 
REAL HA,ANG,ZMAX,DZE 
INTEGER NP 
PI=3.1415926 
BK=2.*PI 
WRITE(6,*) ENTER THE CONE HALF ANGLE IN DEGREES’ 


93 


READ(5,*)HA 
IF(HA.GT.90.)THEN 
WRITE(6,*)’RANGE OF ANGLE: O< ANGLE<90’ 
WRITE(6,*) RANGE EXCEEDED’ 
GOTO 10 
ENDIF 
ANG=HA*PI/180. 
WRITE(6,*) ENTER THE MAXIMUM Z’ 
READ(5,*) ZMAX 
WRITE(6,*)’ ENTER THE NUMBER OF POINTS (NP)? 
READ(S,*)NP 
ZH(1)=ZMAX*BK 
RH(1)=0. 
DZH=ZH(1)/FLOAT(NP-1. ) 
DO 10 I=1,NP 
IF(ZH(I).EQ.0.) THEN 
ZH(I)=.00000001 
ENDIF 
ZH(I+1)=ZH(I)-DZH 
RH(I+1)=RH(1)+(DZH*SIN(CANG) ) 
10 CONTINUE 
BASE=RH (NP) /BK 
ZP=0. 
RS=0. 
RETURN 
END 
CHEER KEKKKKEKKKEKEEKEKKKKEKKKKEKKEKEKEKEXKEEKEEKSUBROUTINE DISK 
SUBROUTINE DISK(NP,ZH,RH,BASE,RS,ZP) 
REAL ZH(400) ,RH(400) ,BASE,RS,ZP,RMAX 
INTEGER NP 
WRITE(6,*)’ENTER THE RADIUS OF THE DISK’ 
READ(5,*)RMAX 
WRITE(6,*)’ENTER THE DISTANCE TO THE DISK’ 
READ(5,*)ZP 
WRITE(6,*)’ENTER THE NUMBER OF POINTS (NP)’ 
READ(5,*)NP 
PI=3.1415926 
BK=2.*PI 
DRH=RMAX*BK/FLOAT(NP~-1. ) 
DO 10 I=1,NP 
ZH(1)=ZP*BK 
RH(I)=DRH*FLOAT(I-1. ) 
10 CONTINUE 
BASE=RMAX 
RS=0. 
RETURN 
END 
CHER EKKEKKEEEEEEEERKEEEEEEEEEEEEEKEKEKEKEKEKEKESUBROUTINE PARAB 
SUBROUTINE PARAB(NP,ZH,RH,DM,FOD,ZP) 
DIMENSION RH(400) ,ZH(400) 
WRITE(6,*) ’ENTER DIAMETER AND F/D RATIO’ 
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READ(S,*) DM,FOD 
WRITE(6,*) ’FEED DISTANCE FROM FOCUS (+ IS NEARER)’ 
READ(S5,*) ZP 
WRITE(6,*) °ENTER NUMBER OF GENERATING POINTS’ 
READ(5,*) NP 
BK=2.*3.14159 
FM=FOD*DM 
PHIV=2.*ATAN(1./4./FOD) 
DO 24 I=1,NP 
TH=FLOAT(I-1)*PHIV/FLOAT(NP-1) 
RM=2.*FM/(1.+COS(TH) ) *BK 
ZH(I)=RM*COS(TH) +ZP*BK 
RH(I)=RM*SIN(TH) 
24 CONTINUE 
RETURN 
END 
CHEEK KEK KKK EEK KKK KEKE KEKEKEKEK KE KEXKHKEKEKEUNCTION TAPER 
FUNCTION TAPER( RHO) 
C SPECIFY ANTENNA ILLUMINATION FUNCTION. REAL FUNCTION OF 
C RHO ONLY 
TAPER=1. 
RETURN 
END 
CHEEK KKK KEK KKK EEE K KEKE KE KEKE KEKE EKEKEKKES*EKS[BROUTINE ANTFF 
SUBROUTINE ANTFF(THETA,PHI,THETAS,PHIS,ARAD,ETF,EPF) 
C COMPUTES CLOSED FORM ANTENNA PATTERN IN THE FAR FIELD. 
C ETF AND EPF ARE RETURNED. ANGLES ARE IN RADIANS. 
C PRESENT ENTRY IS FOR: 
C 
C UNIFORM ILLUMINATION SCANNED TO A DIRECTION (THETAS,PHIS) 
Cc 
COMPLEX ETF,EPF,JK,SCL,EB 
JK=(0.0,6.283185307) 
BK=6 . 283185307 
PI=BK/2. 
C SET RADIATION PATTERN TO ZERO IN THE REAR HEMISPHERE 
IF(COS(THETA).LT.O.) THEN 
EB=(0..,0.) 
ELSE 
BB=SIN( THETA) *SIN(PHI)-SIN(THETAS) *SIN(PHIS) 
AA=SIN(THETA) *COS(PHI)-SIN(THETAS) *COS(PHIS) 
ARG=BK*ARAD*SQRT (AA**2+BB**2) 
SCL=-(0.,1.)*2.*PI**2 
IF(ABS(ARG) .LT.1.E-5) THEN 
EB=SCL/2. 
ELSE 
EB=SCL*(BESSJ1 (ARG) /ARG) 
ENDIF 
ENDIF 
C EB IS THE X COMPONENT; GET ETHETA (ETF) AND EPHI (EPF) 
C ASSUME AN ELEMENT FACTOR, ELMT 


ELMT=SQRT( ABS(COS(THETA) ) ) 
ELMT=1. 
ETF=EB*COS (THETA) *COS(PHI) *ELMT 
EPF=-EB*SIN(PHI)*ELMT 
RETURN 
END 
FUNCTION BESSJ(N,X) 
RETURNS THE BESSEL FUNCTION B OF ORDER N (>1) AND REAL 
ARGUMENT X. 
PARAMETER (IACC=40, BIGNO=1.E10,BIGNI=1.E-10) 
IF(N.EQ.0) THEN 
IF N=0 CALL BESSJO 
BESS J=BESSJO(X) 
ELSE IF(N.EQ.1) THEN 
IF N=1 CALL BESSJ1 
BESS J=BESSJ1(X) 
ELSE 
IF N>1 USE RECURSION 
BESS J=0. 
IF(X.NE.O.) THEN 
TOK= 26 x 
IF(X.GT.FLOAT(N)) THEN 
BJM=BESSJO(X) 
BJ=BESSJ1(X) 
DO 11 J=1,N-1 
BJP=J*TOX*BJ-BJM 
BJIM=BJ 
BJ=BJP 
11 CONTINUE 
BESS J=BJ 
ELSE 
M=2*((N+INT(SQRT(FLOAT(IACC*N) )))/2) 
BESSJ=0. 
JSUM=0. 
SUM=0. 
BJP=0. 
BJ=1. 
DO 1253=M 1-1 
BJM=J*TOX*BJ-BJP 
BJP=BJ 
BJ=BJM 
IF(ABS(BJ).GT.BIGNO) THEN 
BJ=BJ*BIGNI 
BJP=BJP*BIGNI 
BESS J=BESSJ*BIGNI 
SUM=SUM*BIGNI 
ENDIF 
IF(JSUM.NE.0) SUM=SUM+BJ 
JSUM=1-JSUM 
IF(J.EQ.N) BESSJ=BJP 
12 CONTINUE 
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SUM=2. *SUM-BJ 

BESS J=BESSJ/SUM 

ENDIF 

ENDIF 

ENDIF 

RETURN 

END 

FUNCTION BESSJO(X) 
Cc 
C BESSEL FUNCTION OF O ORDER, REAL ARGUMENT X 
C (SEE ’NUMERICAL RECIPES’, P.172) 


c 
REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,Ri,R2,R3,R4,R5,R6, 
* $1,82,S3,54,S5,S6 
DATA P1,P2,P3,P4,P5/1.D0,-.109862827D-2, .2734510407D-4, 
* —.2073370639D-5, .2093887211D-6/ 
DATA Q1,Q2,Q3,04,Q5/-.1562499995D-1, .1430488765D-3, 
* —,6911147651D-5, .7621095161D-6,-.934945152D-7/ 
DATA R1,R2,R3,R4,R5,R6/57568490574. DO ,-13362590354.D0O, 
* 651619640.7D0,-11214424. 18D0, 77392.33017D0, -184.9052456D0/ 
DATA S1,S2,S3,S4,55,S6/57568490411.D0, 1029532985 .D0, 
* 9494680.718D0,59272.64853D0, 267.8532712D0,1.D0/ 
IF(X.EQ.0.) THEN 
BESSJO=1. 
ELSE IF(ABS(X).LT.8.) THEN 
Y=X**2 
BESS JO=(R1+Y*(R2+Y*(R3+Y* (R4+Y*(R5S+Y*R6)))))/ 
* (S1+Y¥*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) 
ELSE : 
AX=ABS(X) 
Z=8./AX 
Y=Z**2 
XX=AX-.785398164 
BESS JO=SQRT(.636619772/AX) *(COS(XX) *(P1+Y*(P2+Y* (P3+ 
* Y*(P4+Y*P5))))-Z*SINCXX) *(Q1+Y*(Q2+Y*(Q3+ 
* Y¥*(Q4+Y*Q5))))) 
ENDIF 
RETURN 
END 
FUNCTION BESSJ1(X) 
€ 


C BESSEL FUNCTION B OF ORDER 1, REAL ARGUMENT X 
C (SEE ’NUMERICAL RECIPES’, P.173) 
Cc 
REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 
* $1,52,53,84,55,S6 
DATA P1,P2,P3,P4,P5/1.D0, .183105D-2,-.3516396496D-4, 
* .2457520174D-5,-.20337019D-6/ 
DATA Q1,Q2,Q3,Q4,Q5/.04687499995D0, - .2002690873D-3, 
* .8449199096D-5,-.99228987D-6, . 105787412D-6/ 
DATA R1,R2,R3,R4,R5, R6/72362614232.D0,-7895059235.D0, 


om 


* 242396853.1D0,-2972611.439D0, 15704. 4826D0, -30. 16036606D0/ 
DATA S1,S2,S3,S4,S5,S6/144725228442 .DO, 2300535178.D0, 
* 18583304. 74D0, 99447 .43394D0 ,376.9991397D0,1.D0/ 
IF(X.EQ.0.) THEN 
BESSJ1=0. 
ELSE IF(ABS(X).LT.8.) THEN 
Y=X*"*2 
BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(RS+Y*RE)))))/ 
* (S1+Y¥*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ))) 
EESE 
AX=ABS(X) 
Z=8./ 8X 
Y=Z**2 
XX=AX-2.356194491 
BESSJ1=SQRT(.636619772/AX)*(COS (XX) *(P1+Y*(P2+Y*(P3+ 
* Y*(P4+Y*P5))))-Z*SIN (XX) *(Q1+Y*(Q2+Y*(Q3+ 
* Y¥*(Q4+Y*Q5)))))*SIGN(1. ,X) 
ENDIF 
RETURN 
END 
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APPENDIX C. GAIN.F LISTING 


Material on the following pages 1s a listing of the program GAIN.F described 


in chapter II]. 


oe 


Cee i RR hk OO EEE EE EEE 


PROGRAM [aint 
DATE : 17 November 1992 
PROGRAMMERS : Dr. D. Jenn, LT K. Klopp 


Computes the gain of a circular aperature antenna radiating through 

a radome consisting of a body of revolution, by numerically 
integrating electric field surrounding system. E-Field is generated by 
currents distributed on a surface of the body of revolution. The 
currents, and their locations were previously determined using the 
program “radome:£: , and Stored in the file 27 (curcoctch onsu, 


INPUT < f£77curcoefsA/B 


COO a aaa a Oo AQABA aA Aa. G 


CHEEK KEKE EEE EEE KEE EEE EKER EEE EEE EEE EEE EEE EEK EEE EES 


CHARACTER ITH 
CHARACTER*12 CC 
CHARACTER*10 GO 
CHARACTER*8 GNT,GNPHI,CGP 
CHARACTER*14 TPTS,PPTS,PHIPTS 
COMPLEX R(1600) ,LEADTERM,ETT,EPT 
COMPLEX EP,ET,ETF,EPF,ETCOMB,EPCOMB 
COMPLEX EXP1,EXP2,CONJG,CEXP,CMPLX ,CT(30000) , IMP, JK 
COMPLEX U,UC,ET1,EP1,ET2,EP2,ZL0(400) 
G DIMENSION RH(400) ,ZH(400) , IPS(800) 
DIMENSION RH(400) ,ZH(400) 
DIMENSION XT(4) ,AT(4) ,A( 100) ,X( 100) ,XP(100) , AP(100) 
DIMENSION Q(30) ,$(30) 
DIMENSION THRR(1000,1000) , PHRR( 1000) 
C DIMENSION ETSCAT(500) ,EPSCAT(500) 
C DIMENSION EXP(500) , ANG(500) , ECP(500) 
INTEGER CNPHI,SELECTION 
DATA PI/3.1415926/ 
DATA RR/1000/ 
Rad=PI/180. 
BK=2.*PI 
V=(On 
U0=(0. ,0.) 
UC=-U/4./PI 
JK=(0.0,6.283185307) 
C*****ENTER A LETTER TO INDICATE WHICH ITERATION THIS [5**4**«*¢4e4e444%%%% 
Cee + + (ALL DATA FILES ARE APPENDED WITH THIS LETTER) eK EEK KEKE 
Cee * (.m FILES ARE FOR GENERATING PLOTS IN MATLAB) ee 
WRITE(6,*)’ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS’ 
READ(S,*)ITH 
€C="curcoetsdat 7/1 0H 
GO=’gainout’//ITH//’.m’ 
CKEKKKKKKEKEEKE KEKE KEKE KEKE EEE EEE KEKEKKKEEEKKEKEKEEKKEKEKEEEEEEEEEEE KEES 
c 
C READ in the following data from file called f77curcoefs: 


100 


C 

C 1teration, ITH 

C geometry of BOR, SELECTION 
& radius of base of BOR, BASE 

© total points, NP 

e number of t-intervals, MT 

C phi-intervals, ME 

C sum of t and phi intervals, N 

c number of rows of currents, NROW 

C total number of modes, MODES 

C number of sub vectors, NBLOCK 

C antenna scan angle in theta, THS 

C antenna scan angle in phi, PHS 

c antenna radius, ARAD 

C observation angle, PHIO 

C complex radome impedance, 16 Us 

C Zprime, ZP 

C surface radius, RS 

C closed form antenna pattern, IFF=0 

C 

C files where Gauss integration weights and abscissas are stored 
c 

c coordinates of surface points *2PI/lamda, ZH(I), and RH(I) 
c 


CHER KKKEKKE KKK KKKEKREKREKK KK KKK EKKEKKK KEK KEKE EEE EK EEE KK EEE KEKE KKK EEE EEE EKK 
OPEN(1,file=CC) 
WRITE(6,*) ’FILE OPENED IS ’,CC 
READ(1,*) ITH,SELECTION,BASE,NP,MT,MP,N 
READ(1,*) NROW,MODES ,NBLOCK,THS,PHS,ARAD, PHIO 
READ(1,*) IMP,ZP,RS,IFF 


C write(6,*)IMP,ZP,RS,IFF 
READ(1,*) (RE(L),L=1,NP) 

C write(6,*) (RH(L),L=1,NP) 
READ(1,*) (ZH(L) ,L=1,NP) 

Cc write(6,*) (ZH(L),L=1,NP) 
READ(1,*) (ZLO(CL) ,L=1,NP) 

Cc write(6,*) (ZLO(L) ,L=1,NP) 


TTL E LE LESS ESE SSL ESET TESTES SSS TET TEES TC ESE SESE TEESE STC ESET SCTE STE SSCS ETS SES ETS FS 
G 
C READ in current coefficients, CT(I) from file curcoefsdatA/B 
c (This is one one collumn vector.) 
& 
CHEER KEKE KEE KKK KEKE KE EEK KK KEKE EEE EEE EEE EEE EEE EEK KEE EEE ERE KEKE EEE EEEEEEES 
READ(1,*) (CT(L) ,L=1,NROW) 
C write(6,*)(CT(L),L=1,NROW) 
READ(1,*) GNT,GNPHI,CGP 
WRITE(6,*) ’RUN CODE READ FROM DISC : ’,ITH 
THETAS=THS*RAD 
PHIS=PHS*RAD 
MHI=MODES+1 
CLOSE(1) 
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Cte RR RR RRR ORR 


Antenna integration in x,y direction: PHIPTS = gaus/gaus#; XP(I),AP(I) 


C 

C Read in weights A(I), and Abscissas X(I) for the following Gauss- 

C Quadrature Integrations: 

C fa Ke wts absc 
c 

G BOR integration in t direction: TPTS = gaus/gaus#; XT(I),AT(I) 
C BOR integration in phi direction: NPHI = gaus/gaus#; X (1I),A (I) 
C 

G 


STS SSS SSS ESTES ESET ES SCS ESSE SCTE SESS ESTES EEE TESST TCE ETT TE SCS TTS T ST TCS SS ST STS STS SS SS eS 
TPTS=’gaus/’//GNT 
PPTS=’ gaus/’//GNPHI 
PHIPTS=’gaus/’//CGP 
OPEN (2,FILE=TPTS ,STATUS=’ OLD’ ) 
OPEN(3,FILE=PPTS ,STATUS=’OLD’ ) 
OPEN(8,FILE=PHIPTS , STATUS=’OLD’ ) 
READ(2,*)NT 
WRITE(6,*)’NT = ’,NT 
IF(NT.GT.4)THEN 
WRITE(6,*)’?MAXIMUM NUMBER OF POINTS(NT) IS 4? 
GOTO 999 
ENDIF 
READ(3,*)NPHI 
WRITE(6,*)’NPHI = ’,NPHI 
IF(NPHI.GT.200) THEN 
WRITE(6,*) ’MAXIMUM NUMBER OF POINTS(NPHI) IS 200’ 
GOTO 999 
ENDIF 
READ(8,*) CNPHI 
WRITE(6,*)’CNPHI = ’,CNPHI 
IF (CNPHI.GT.100)THEN 
WRITE(6,*)’?MAXIMUM NUMBER OF POINTS(CNPHI) IS 100’ 
GOTO 999 
ENDIF 
C LOAD THE WEIGHTS AND ABSCISSAS IN THE VECTORS. 
DO 1 M=1,NT 
READ(2,*,END=1)XT(M) , AT(M) 
WRITE(6,*)XT(M) ,AT(M) 
1 CONTINUE 
DO 2 M=1,NPHI 
READ(3,*,END=2)X(M) , ACM) 
WRITE(6,*)X(M) ,A(M) 
2 CONTINUE 
DO 4 M=1,CNPHI 
READ (8,*,END=4)XP(M) , AP(M) 
WRITE(6,*)XP(M) ,AP(M) 
4 CONTINUE 
CLOSE(2) 
CLOSE(3) 
CLOSE(8) 


Ce OR OO RO OR OO RK KEE EERE EEE KEE EK REE EEE KE EK KEKE 
g 
C Create an output file called ’gainoutA/B.m’,and output data on geometry etc. 
C 


COR OOOO OR ORO ORO OR OR EEK OK KK 


OPEN(7 , FILE=GO0) 


WRITE(7,2000)’ GEOMETRY = ’,SELECTION,’ NP = ’,NP 
&, (oa? Sl Me Se MPN) = IN 
&, > MODES = ’,MODES,’ NBLOCK = ’ ,NBLOCK 
a ce 0 Ceram) 0 
WRITE(7,2001)’ BASE = ’,BASE,’ ARAD = ’,ARAD,’ THETASCAN = ’,THS 
& > PHISCAN = ’,PHS 


WRITE(7,*)’NPHI =’ ,NPHI 
2000 FORMAT(//5X,’%% RADOME GAIN CALCULATION %%’//,10(/°% °,A13,14)) 
2001 FORMAT(3(/’% °,A13,F9.2),/) 


CHEEK KEKE KEKE EK REE KEKE KEKE KEK KER KEKE KEK KEE KEKE KEKE KEKE KKK EKEKEKEKKEE KE 


Gain Computation: 
4 pi * Max Radiation Intensity 


Integration Over Closed Surface ( Radiation 
Intensity * Beam Solid Angle ) 


Radiation Intensity = (1/(2*eta))*| E(t,phi,r) |**2 
Since E is in the far field, need to only determine Etheta, Ephi. 
|E|**2 = JE scattered|**2 + |E antennal **2 


|E scattered| = Etheta + Ephi 


Gir +O) Git Co Cl CieGi Ol Ca C2 Ca Cl 


CEE EKEKE KEE KKEREKE EKER EKER KEKE KEK KEE KEE EEE EEE EEE EKER EKER EEE EERE EEE EEE KEE KEK 
C*x*Enter # divisions for integration in theta, and phi 
C directions. (This allows easy variation of the number of integration steps 
C without having to recompute a different set of weights and abscissas for 
C Gauss Integration.) 
WRITE(6,*)’ENTER NUMBER OF DIVISIONS IN THETA DIRECTION’, 


& ae Dans 
READ(S,*)NDIVT 
WRITE(6,*)’NDIVT = ’,NDIVT 
WRITE(7,*)’% NDIVT = °,NDIVT 


WRITE(6,*)’ENTER NUMBER OF DIVISIONS IN PHI DIRECTION, NDIVP = 7’ 
READ(5,*)NDIVP 


WRITE(6,*)’NDIVP = ’,NDIVP 
WRITE(7,*)’% NDIVP = ’,NDIVP 
KK=NDIVT 

C do 500 NDIVT=1,kk 


NDIVP=NDIVT 
CeO GOO Ook aoocrraeceintegration intervals 0<theta<180, 0<ph1<360 
THETA1=0 
THETA2=180 
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S1=THETA1 
S2=THETA2 
$1=S1*Rad 
S$2=S2*Rad 
PHI1=0 
PHI2=360 
Q1i=PHI1 
Q2=PHI2 
Q1=Q1*Rad 
Q2=Q2*Rad 
CHEeKEKKKEKEKKEKKEKAREKRARK KEK AKL ARERR ELLR#CHeta integration angles 
DS=(S2-S1)/FLOAT(NDIVT) 
DO 200 I=1,NDIVT+1 
S(T)—=C1=1)*Ds 
200 CONTINUE 
CHHREEEREREKHERE RRA HRERK AR ERHRKKEEREKDHI] Integration angles 
DQ=(Q2-Q1)/FLOAT(NDIVP) 
DO 201 I=1,NDIVP+1 
Q(I)=(I-1)*DQ 
201 CONTINUE 
SUM=0. 
SUMNONE=0O. 
UMAX=0. 
UMAXNONE=O. 
EMAX=0. 
EMAXNONE=0. 
CeeeeeeeeeeeeeeeeeeeeeeeeV phi integration of rad intensity * beam solid angle V 
DO 300 JJ=1,NDIVP 
P1=DQ/2. 
P2=(Q(JJ+1)+Q(JJ))/2. 
CHEEREEEER EERE ER ERE EKA KEREEHE*EY Phi integration determining Escattered V 
C V within each interval of phi, NDIVP V 
DO 300 J=1,NPHI 
PHR=P1*X(J)+P2 
Ceaeeekektekeeeeee)V theta integration of rad intensity * beam solid angle V 
DO 301 II=1,NDIVT 
T1i=DS/2. 
T2=(S(lEtt) +5 (le 2. 
CeeeeaeceeecekeeeeeeceeeceeeyY theta integration determining Escattered V 
C V within each interval of theta, NDIVT V 
DO 301 I=1,NPHI 
THR=T1*X(1I)+T2 
TP(THR. or Fl Gd tO 302 
THR=(BK-THR) 
PHR=PHR+PI 
80Z CONTINUE 
CHEEK EEE EERE KEKE KEKE KERAEKEREEKEREKE RK KKK KKK RKREKBe aM Solid Angle 
STRA=SIN(THR)*A(I) *A(J) 
KJ=(JJ-1)*NPHI+J 
KI=(II-1)*NPHI+I 
THRR(KJ ,KI)=THR*180/PI 
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PHRR(KJ)=PHR*180/PI 
SNP=SIN(PHR) 
CNP=COS(PHR) 


CR KEKE KEE KEKE EEE EERE EERE EE EEE EERE EEE EEKKKE EEE 


C _ ne = = = s 

oy | | I | 

c | s | [ te theta Sphistueta |) =-t> -| 
SE | -jkr MHI | R R ie ae | 

cro tl -jhe aeo | n | | on  [ jnphi 
S| =) . eee el lea le 

eo | s | 4pir j= | Stearns = pad ilo Ohi 

Se ie Cl n=1 | R R (es | 

Seyi phil a n (ote 
Sewia. ie _l ote =| 


CU EREEREEEEEREREKESEEEAEKERERKRERESCG Page 25. H,E, and Combined Field*«*«x*x«V 
© -jkr 
C —jhe 


C Leading Term eta (h) is cancelled by its inverse in CT (I), 


C eJophg r dependence is ignored, and the whole Term is 
C scaled by 2pi for gauss quadrature integration 
C 


C Leading Term is reduced to (-j/4pi)/2pi 
Co OO RII III III III EOI ROY 


c LEADTERM=UC/BK 
LEADTERM=UC 
ET=(0,0) 
EP=(0,0) 
ET1=(0.,0.) 
EP1=(0. ,0.) 
Bo (OneOr) 
EP2=(0.,0.) 

CREEK K KEKE KKK KKKKEKKKKEKKKEKKEKEKKKEEK KEKE EK KKK KKKKKKKKK KEKE KEKE KEKE KEKE EK! 

c MHI 

C I 

C \ 

Cc jo 

C n=1 


CHEEK KK KKKKKKKKKKKKKKKKEKEKEKEKKEKEKEKEKKKEEKKEKKEKKKEKE KEKE KKK KEK EKKEKEKEKKKEKEEEEKEEEY 
DO 305 M=1,MHI 
NM=M-1 
CHEER E EEE EERE EERE EERE EEE EE EE EEE EEE EEE EERE EKER EERE KER EKER EEA JN DHI 
C e 
EXP1=CEXP(CMPLX(0O. , FLOAT(NM) *PHR) ) 
EXP2=CONJG(EXP1) 
Cxeaeeeeeeeeeedetermine plane wave excitation vector in far field, R 
CALL PLANE(NM,NM,NP,NT,RH,ZH,XT,AT,THR,R) 
NTOP1=MODES-NM 
NTOP 2=NBLOCK-(NTOP1+1) 
NS2=NTOP1*N 
NS1=NTOP2*N 


CHEEK KKKKKKKKKKKKKKKKKKKEKKEKKEKKEKKEEKEKKKKKEKKKEKKEKEEKKEKKKEKEKKEKKKEEKEE 
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¢ | . t,thetal | _t | jmphi | 2 t) phi Jie Gens 
€ iE: heap | e vand || R (ere aire 
c | on [oil Waeiriiass AI [om |) ace al 
CREEK EEKEREKEKREKEKEEEREREREREREEKEKEKERERERERERERERERERERERERKEK EK 
DO 306 L=1,MT 
CREE EEE EE ERE EERE EEE EEE EERE EEE EEREREKEEKDOSI tI VE mode +n 
ET1=ET1+R(L) *CT(L+NS1) *EXP1 
EP1=EP1+R(L+N) *CT(L+NS1) *EXP1 
CHEEK ERE REE EEE EERE EEA EA EAHA EE EKEEKEEEKESEN OTATIVE mode -n 
IF(NM.EQ.0) GO TO 306 
ET1=ET1+R(L) *CT(L+NS2) *EXP2 
EP1=EP1-R(L+N) *CT(L+NS2) *EXP2 
306 CONTINUE 


C¥EEEKKEXEXKEXKKEKEKEKEKEKKKEKKKEEKEKEKERKEEKEKEAREKAERERKEAKEEEKEREEEKE KE 

C lepha Cheta | ale to ai npha |) Se pha paves el a | egnens 

Cie ee ee l e ,and | R lee | e 

© ilern [Sere | (ren a: 

C¥EEREKKXEKKXEKXKXEXKXEKXKEKEKEEREKEKEKEREKEKEKKREKRKKREKEKEKREKKKKEKEEEKEKEKE 

DO 307 L=1,MP 

CHREREXEERERERAEREERERE LAE RAEAEKAEKEEEREREELAKEADOSI CLIVE mode +n 
ET2=ET2+R(L+MT) *CT(L+NS1+MT) *EXP1 
EP2=EP2+R(L+MT+N)*CT(L+NS1+MT ) *EXP1 

C¥SKAEEEEERSESES EELS“ SEASELELE SESSA SET ESRESEES ESE ET COAT IVE mode -n 
IF(NM.EQ.0) GO TO 307 
ET2=ET2-R(L+MT) *CT (L+NS2+MT) *EXP2 
EP2=EP2+R(L+MT+N)*CT(L+NS2+MT ) *EXP2 


307 CONTINUE 

é ET=ET+LEADTERM*(ET1+ET2) 
G EP=EP+LEADTERM*(EP1+EP2) 
305 CONTINUE 


Cok EKER EEKEEKEKEKK™ EFC mode summation ~ 
ETT=LEADTERM*(ET1+ET2) 
EPT=LEADTERM*(EP1+EP2) 
ea KEKE EEKEKEKEKEKEKKEKKEEKOGG in antenna contribution 
C FEED CONTRIBUTION IN THE FAR FIELD IS ETF,EPF 
Ca a REE EEE EKEEKEXKB_ESS J 1. 
C USE CLOSED FORM FEED EXPRESSION FROM SUBROUTINE ANTFF IF 
C IFF=0; OTHERWISE USE BRUTE FORCE EVALUATION FROM CIRCRTP 
C 
ROBS=1000. 
IF(IFF.EQ.0) THEN 
CALL ANTFF(THR,PHR,THS,PHS,ARAD, ETF, EPF) 
€ write(6,*)” ET = ° CABS(EIy ETF = "(CABS @EIE) 
ELSE 
RHB=ROBS*SIN(THR) 
ZHB=ROBS*COS (THR) 
CALL CIRCRTP(CNPHI,XP,AP,ARAD,THS, PHS, 
& PHR ,RHB,ZHB,CIRCR,CIRCT,CIRCP) 
C REMOVE THE 1/R DEPENDENCE BECAUSE EXP(-jKR)/R IS OMITTED IN 
C THE SCATTERED FIELDS ET AND EP 
ETF=CIRCT*ROBS 
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EPF=CIRCP*ROBS 


Cc write(6,*)’ ET = ’,CABS(ET),’ ETF(ROBS) = ’,CABS(ETF) 

ENDIF 
C¥ERKKREEKEKEKERERRAEKERAEEARERE EERE AAEAEKEAKESREATOCAL E theta and E phi 
C auitveco,*)) THR(’))),Js21,2,NM,°) = °,180*THR/pi 


ETCOMB=ETT+ETF 
EPCOMB=EPT+EPF 
EMAG=SQRT (CABS (ETCOMB )**2+CABS (EPCOMB) **2) 
EMAGNONE=SQRT(CABS (ETF) **2+CABS (EPF) **2) 
IF (EMAG.GT.EMAX) THEN 
EMAX=EMAG 
EMAXNONE=EMAGNONE 
WRITE(6,*)’ EMAG(’,J,I,’) = ’,EMAG 
END IF 
CHEEKREKEEEEEEEEREEKAEEEEAEKAKKAREKEKAREEKEKEAGCETerMine max field intensity 
UMAG=EMAG**2 
UMAGNONE=EMAGNONE**2 
IF (UMAG.GT.UMAX ) THEN 
UMAX=UMAG 
THETAMAX=180*THR/PI 
PHIMAX=180*PHR/PI 


C WRITE(6,*)’ UMAX = ’,UMAX 
© WRITE(6,*)’ THETAMAX = ’,THETAMAX 
Cc WRITE(6,*)’ PHIMAX = ’, PHIMAX 

END IF 


IF (UMAGNONE. GT . UMAXNONE ) THEN 
UMAXNONE=UMAGNONE 
ENDeIy 
CeeeeeeeeeeeeeeeeRARERREMULtIply integration. by beam solid angle (STRA) 
SUM=SUM+UMAG*STRA 
SUMNONE=SUMNONE+UMAGNONE*STRA 
301 CONTINUE 
CEE KAEAREKKAKEEKE™ CN theta integration - 
300 CONTINUE 
CHEEK AEEREEEKERKEAEEKEKREREEKAEEEKEKKAEKEREEEAEEEEEKE” CNC phi integration es 
CHEREEEEKKKKERERRKREKEREERESCale each summation by (b-a)/2 for Gauss integration 
PRAD=T1*P1*SUM 
PRADNONE=T1*P1*SUMNONE 
WRITE(6,*)’PRAD = ’,PRAD 
WRITE(6,*)’PRAD NO RADOME = ’,PRADNONE 
OCSTSCESESSESESE SESS SESE SS SSS SS SSS SS SS SESS STS SSS SS SST S SS SSS SS TT ESE STS SESS ST ES SS SS 
C End Integration over closed surface radiation intensity * beam solid angle 
CFE SESES ESE SESE SES SEES EES TSS SESE SESE SESS ESTEE SELES SESE SES SE SES EES SE TSS SET ES ES ES SS 
CHEEK ERRE RARER EERE REE KE RAEEA EKER RARER RRA ERRAERAEAARA AAA RAAKAKCOMPUTE Galn 
GAIN=4*PI*UMAX/PRAD 
GAINNONE=4*PI*UMAXNONE/PRADNONE 
GDB=10. *alog10(GAIN) 
GDBNONE=10. *alog10(GAINNONE) 
WRITE(6,*)’GAIN’ ,NDIVT,’ = ’,GAIN,’ IN DB = ’,GDB 
WRITE(6,*)’GAINNONE = ’,GAINNONE,’ IN DB = ’,GDBNONE 
WRITE(7,*)’NDIVT(’,NDIVT,’)=’,NDIVT,’; NDIVP(’ ,NDIVP,’)=’,NDIVP 


Oe 


WRITE(7 ,2003)’GAIN(’ ,NDIVT,’ )=’ ,GAIN,’;GAINDB(’ ,NDIVT,’ )=’ ,GDB 
WRITE(7 ,2003)’GAINNONE(’,NDIVT,’) = ’,GAINNONE, 
&’;:GAINNONEDB(’ ,NDIVT,’) = ’,GDBNONE 
2003 FORMAT(/(A, 1e¢hs Fb te Seis nF 1025) 
C500 continue 
STOP 
999 end 
Ce ee EEE EE EEK EKER KEKE KE REKKKKKKEEEKKE*KEE ND OF MAIN PROGRAM 
Ce EK KKK EEE KEKE KKKKKEKKEKKKKKKKKEEKEEXSUBROUTINE PLANE 
C REFERENCE: TECHNICAL REPORT TR-80-1 ;(pages 57-64) 
SUBROUTINE PLANE(M1,M2,NP,NT,RH,ZH,XT,AT,THR,R) 


PLANE WAVE EXCITATION VECTOR IN THE FAR FIELD. FROM MAUTZ AND HARRINGTON. 
MODIFIED TO DO ONLY ONE ANGLE AND FREQUENCY PER CALL. 
*x* NEW BESSEL FUNCTION ROUTINE FROM NUM. RECIPES ** 


A OO) OO) 


COMPLEX R(1600) ,U,U1,UA,UB,FA(1500) ,FB(1500) ,F2A,F2B,F1A,F1iB 
COMPLEX U2,U3,U4,US5 
DIMENSION RH(400) ,ZH( 400) ,XT(4) , AT(4) ,R2(4) ,Z2(4) 
MP=NP-1 
MT=MP-1 
N=MT+MP 
N2=2*N 
CC=COS (THR) 
SS=SIN (THR) 
U=(O a) 
U1=3.141593*U**M1 
M3=Mi+1 
M4=M2+3 
IF(M1.EQ.0) M3=2 
M5=M1+2 
M6=M2+2 
DO 12 IP=1,MP 
K2=IP 
IT=IP+1 
DR=.5*(RH(I)-RH(IP)) 
DZ=.5*(ZH(1)-ZH(IP)) 
D1=SQRT(DR*DR+DZ*DZ) 
Ri=.25*(RH(I)+RH(IP) ) 
IF(R1.EQ.0.) Ri=1. 
Zi=.5*(ZH(1)+ZHCIP)) 
DR=.5*DR 
D2=DR/R1 
DO 13 L=1, Nm 
R2(L)=R1+DR*XT(L) 
Z2(L)=Z1tDZex7T ce) 

13 CONTINUE 
D3=DR*CC 
D4=-DZ*SS 
D5=D1*CC 
DO 23 M=M3,M4 
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FA(M)=0. 
FB(M)=0. 

23 CONTINUE 
BO 15 L=1 NT 
X=SS*R2(L)*2. 
ARG=Z2(L) *CC 
UA=AT(L)*CMPLX(COS(ARG) , SIN( ARG) ) 
UB=XT(L)*UA 

C THIS LINE REPLACES HARRINGTON’S BLOCK 

DO 25 M=M3,M4 
BES=BESSJ(M-2,X) 
FA(M)=BES*UA+FA(M) 
FB(M)=BES*UB+FB(M) 

25 CONTINUE 

15 CONTINUE 
IF(M1.NE.0) GO TO 26 
FA(1)=-FA(3) 
FB(1)=-FB(3) 

26 UA=U1 
DO 27 M=MS5,M6 
M7=M-1 
M8=M+1 
F2A=UA*(FA(M8)+FA(M7) ) 
F2B=UA*(FB(M8)+FB(M7) ) 
UB=U*UA 
F1A=UB*(FA(M8)-FA(M7) ) 
F1B=UB*(FB(M8)-FB(M7) ) 
U4=D4*UA 
U2=D3*F1A+U4*FA(M) 
U3=D3*F1B+U4*FB(M) 
U4=DR*F2A 
US=DR*F2B 
K1=K2-1 
K4=Ki+N 
KS=K2+N 
R(K2+MT)=-D5*(F2A+D2*F2B ) 
R(K5+MT)=D1*(F1A+D2*F1B) 
Fe(re-EG-1) GO TO 21 
R(K1)=R(K1)+U2-U3 
R(K4)=R(K4)+U4-U5 
IF(IP.EQ.MP) GO TO 22 

21 R(K2)=U2+U3 
R(KS)=U4+US 

22 K2=K2+N2 
UA=UB 

27 CONTINUE 

12 CONTINUE 
RETURN 
END 

C¥EEEKEKEXKEEKEEEKEXEKEKEKEKEEKEKESKEEEEKEKE*XEXESUITRBROUTINE ANTEFF 

SUBROUTINE ANTFF(THETA,PHI,THETAS,PHIS,ARAD,ETF,EPF) 
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COMPUTES CLOSED FORM ANTENNA PATTERN IN THE FAR FIELD. 
ETF AND EPF ARE RETURNED. ANGLES ARE IN RADIANS. 
PRESENT ENTRY IS FOR: 


UNIFORM ILLUMINATION SCANNED TO A DIRECTION (THETAS,PHIS) 


Cy 1G Cr Cy “Ch CG) 


COMPLEX ETF,EPF,JK,SCL,EB 

JK=(0.0,6.283185307) 

BK=6. 283185307 

PI=BK/2. 

C SET RADIATION PATTERN TO ZERO IN THE REAR HEMISPHERE 

IF (COS( THETA) -ET.0. )) THEN 
EB=(0.,0.) 

ELSE 
BB=SIN(THETA)*SIN(PHI)-SIN(THETAS) *SIN(PHIS) 
AA=SIN(THETA)*COS(PHI)-SIN(THETAS) *COS(PHIS) 
ARG=BK *ARAD*SQRT (AA**2+BB**2) 
SCL=-(0.,1.)*2.*PI**2 
IF(ABS(ARG).LT.1.E-5)THEN 

EB=SCL/2. 
ELSE 
EB=SCL*(BESSJ1(ARG)/ARG) 
ENDIF 
ENDIF 
C EB IS THE X COMPONENT; GET ETHETA (ETF) AND EPHI (EPF) 
C ASSUME AN ELEMENT FACTOR, ELMT 
ELMT=SQRT(ABS(COS( THETA) ) ) 
ETF=EB*COS (THETA) *COS (PHI) *ELMT 
EPF=~EB*SIN(PHI) *ELMT 
RETURN 
END 
Ci i i KEE KEE EEK EERE EKKEKEXKXRESSEL FUNCTIONS 
FUNCTION BESSJ(N,X) 
C RETURNS THE BESSEL FUNCTION B OF ORDER N (>1) AND REAL 
C ARGUMENT X. 
PARAMETER (IACC=40,BIGNO=1.E10,BIGNI=1.E-10) 
IF(N.EQ.0) THEN 
C IF N=0 CALL BESSJO 
BESSJ=BESSJO(X) 
ELSE IF(N.EQ.1) THEN 
C IF N=1 GARE BESSJ1 
BESS J=BESSJ1(X) 
ELSE 
C IF N>1 USE RECURSION 

BESSJ=0. 

IF(X.NE.O.) THEN 

TO x 

IF(X.GT.FLOAT(N)) THEN 

BJM=BESSJO(X) 

BJ=BESSJ1(X) 

DO 11 J=1,N-1 
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BJP=J*TOX*BJ-BJM 
BIM=BJ 
BJ=BJP 
11 CONTINUE 
BESSJ=BJ 
ELSE 
M=2* ((N+INT(SQRT(FLOAT(IACC#N) )))/2) 
BESSJ=0. 
JSUM=0. 
SUM=0. 
BIP=0. 
Bi- 1 
DO 12 J=M,1,-1 
BJM=J*TOX*BJ-BJP 
BJP=BJ 
BJ=BJM 
IF(ABS(BJ).GT.BIGNO) THEN 
BJ=BJ*BIGNI 
BJP=BJP*BIGNI 
BESS J=BESSJ*BIGNI 
SUM=SUM*BIGNI 
ENDIF 
IF(JSUM.NE.O) SUM=SUM+BJ 
JSUM=1-JSUM 
IF(J.EQ.N) BESSJ=BJP 
12 CONTINUE 
SUM=2.*SUM-BJ 
BESS J=BESS J/SUM 
ENDIF 
ENDIF 
ENDIF 
RETURN 
END 
FUNCTION BESSJO(X) 
C 
C BESSEL FUNCTION OF O ORDER, REAL ARGUMENT X 
C (SEE ’NUMERICAL RECIPES’, P.172) 
C 
REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 
* $1,52,53,S4,S5,S6 
DATA P1,P2,P3,P4,P5/1.D0,~-.109862827D~-2, .2734510407D-4, 
* - ,2073370639D-5, .2093887211D-6/ 
DATA Q1,Q2,Q3,04,Q5/-.1562499995D-1, .1430488765D-3, 
* -.6911147651D-5, .7621095161D-6,- .934945152D-7/ 
DATA R1,R2,R3,R4,R5,R6/57568490574. D0, -13362590354.D0, 
* 651619640. 7D0,-11214424. 18D0, 77392. 33017D0, -184.9052456D0/ 
DATA S1,S2,S83,S84,S5,S6/57568490411.D0,1029532985.D0, 
* 9494680.718D0 ,59272.64853D0, 267.8532712D0,1.D0/ 
IF(X.EQ.0.) THEN 
BESS JO=1. 
ELSE IF(ABS(X).LT.8.) THEN 


ey 


Y=Xee2 
BESS JO=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(RS+Y*R6)))))/ 
* (S1+¥*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ))) 
ELSE 
AX=ABS(X) 
Z=8./AX 
Y=Z**2 
XX=AX-. 785398164 
BESS JO=SQRT ( .636619772/AX) *(COS(XX) *(P1+Y¥*(P2+Y*(P3+ 
* Y*«(P4+Y*P5) )))-Z*SIN(XX) *(Q1+Y*(Q2+Y*(Q3+ 
* Y*(Q4+Y*Q5))))) 
ENDIF 
RETURN 
END 
FUNCTION BESSJ1(X) 
C 
C BESSEL FUNCTION B OF ORDER 1, REAL ARGUMENT X 
C (SEE ’NUMERICAL RECIPES’, P.173) 
E 
REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 
* S$1,S$2,S3,S4,S5,S6 
DATA P1,P2,P3,P4,P5/1.D0, .183105D-2,-.3516396496D-4, 
* .2457520174D-5, -.20337019D-6/ 
DATA Q1,Q2,03,Q4,Q5/ .04687499995D0, -. 2002690873D-3, 
* .8449199096D-5,-.99228987D-6, . 105787412D-6/ 
DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235 .D0, 
* 242396853 .1D0,-2972611.439D0, 15704. 4826D0 , -30. 16036606D0/ 
DATA S1,S$2,S3,S4,S5,S6/144725228442.D0O, 2300535178.D0, 
* 18583304. 74D0 , 99447 .43394D0 ,376.9991397D0, 1.D0/ 
IF(X.EQ.0.) THEN 
BESSJ1=0. 
ELSE IF(ABS(X).LT.8.) THEN 
Y=X**2 
BESS J1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(RS+Y*R6)))))/ 
* (S1+Y* (S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) )))) 
ELSE 
AX=ABS(X) 
Z=8e 7k 
Y=Z**2 
XX=AX-2.356194491 
BESS J1=SQRT( .636619772/AX) *(COS (XX) *(P1+Y*(P2+Y*(P3+ 
* Y*(P4+Y*P5) )))-Z*#SIN(XX)*(Q1+Y*(Q2+Y*(Q3+ 
* Y¥*(Q4+Y*QS) ))))*SIGN(1.,X) 
ENDIF 
RETURN 
END 


CEERKEEEKKKKKEEREK EKER REE EEE EEE EEEEEEKKEEKEKEKEKEESITRBROUTINE CIRCRTP 
C SUBROUTINE : CIRCRTP 
C DATE : 1 OCTOBER 1991 


eZ 


REVISED 26 February 1992 
PROGRAMMER R.M. FRANCIS 
2ND REVISION : 27 October 1992 
PROGRAMMER : Kk. A. KECPP 


COMMENTS : THIS SUBROUTINE WILL RIGOROUSLY CALCULATE THE ELECTRIC FIELD 
FOR A CIRCULAR APERTURE. THE FIELD IS CALCULATED AT THE COORDINATES 
SPECIFIED BY RH(.) AND ZH(.) THE APERTURE IS LOCATED AT Z = 0, 

AND SCANNED TO A DIRECTION (THS,PHS). THE SUBROUTINE OGIVE IS THE 

SOURCE OF THE GEOMETRIC DATA REQUIRED BY CIRC TO PERFORM COMPUTATIONS. 
ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 
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Coke RRR III OOO OR OO IK OOK OK 
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1 x \ xt 
E = -jn/(4*pi*k) /_ | ar tee ere )*TAPER(rp) 
x=0 y=0 
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CHEEK KKEKKKKKKKEKKEKEKKKEKEKKEKKKE KKK EKER EKER KKK KKK EEKEKEEK KKK KEKE 
SUBROUTINE CIRCRTP(CNPHI,XP,AP,ARAD,THS,PHS, 
* PHI ,RZ,Z,CIRCR,CIRCT ,CIRCP) 
INTEGER CNPHI 
DIMENSION XP(100) ,AP(100) 
REAL K 
COMPLEX J,JK,G1,G2,X1,CON,CC,SUMX,SUMY,SUMZ 
COMPLEX CIRCR,CIRCT,CIRCP 
PI=3.141592654 
ETA=120.*PI 
C SET FIELD COMPONENTS TO ZERO IN THE REAR HEMISPHERE 
CIRCR=(0,.0-) 
CIRCT=(0. ,0.) 
CIRCP=(0. ,0.) 
C Z OF THE OBSERVATION POINT IS >= 0 IN THE FORWARD HEMISPHERE 
IF(Z.GE.0.) THEN 
SUMX=(0.0,0.0) 
SUMY=(0.0,0.0) 
SUMZ=(0.0,0.0) 
K=6. 283185307 
Jc= (0.0, 1-0) 
JK= (0.0,6.283185307) 
CON=-J/(4.*PI) 
C OMIT THE FACTOR OF ETA SINCE IT IS NOT INCLUDED IN SUBROUTINE ZMAT. 
ARAD**2 IS FOR SCALING OF GAUSSIAN INTEGRATION. READJUST SCALING 
C TO AGREE WITH OLD SUBROUTINE SPHERE (ARBITRARY) 
EC=..CON 
C OUTER LOOP: INTEGRATE OVER X, -ARAD < X < ARAD. 
INNER LOOP: INTEGRATE OVER Y, -ARAD < Y < ARAD. 
C EVALUATE ANTENNA FIELD AT RZ,Z: 
XO=RZ*COS (PHI) 
YO=RZ*SIN (PHI) 


QO? 


‘?) 


I 


STS=SIN(-THS) 
S=RZ/SQRT(RZ**2+Z**2) 
C=Z/SQRT(RZ**2+Z**2) 


C INTEGRATE IN X,Y COORDINATES 


DO 50 M=1,CNPHI 
X=ARAD*XP(M) 

DO 60 N=1,CNPHI 
Y=ARAD*XP(N) 
RP=SQRT (X**2+Y*¥*2) 


C SKIP INTEGRATION POINTS BEYOND THE ANTENNA RADIUS 


60 
50 


IF(RP.LE.ARAD) THEN 
PHIPRIME=ATAN2(Y ,X+1.E-10) 
R=SQRT(((RZ-RP) ##2+Z**2)+(4*RZ*RP*SIN( (PHI-PHIPRIME) /2)**2) ) 
Gi=(((K*R) **2)-1-(JK*R) ) /R**3 
G2=(3+(3*IK*R)-(K*R) **2) /R*¥*5 
X1=JK*(RP*COS(PHIPRIME-PHS) *STS-R) 
Y1=X0-X 
X2=Y1i**2 
Y2=YO-Y 
SUMX=SUMX+(CC*AP(M) *AP(N) *(G1+(X2*G2) ) *CEXP(X1) )*TAPER(RP) 
SUMY=SUMY+(CC*AP(M)*AP(N)*Y1*Y2*G2*CEXP(X1))*TAPER(RP) 
SUMZ=SUMZ+(CC*AP(M) *AP(N) *Y1*Z*G2*CEXP(X1))*TAPER (RP) 
ENDIF 
CONTINUE 
CONTINUE 


C ASSUME AN ELEMENT FACTOR, ELMT 


ELMT=SQRT(ABS(C) ) 

CIRCT=(C*COS (PHI )*SUMX+C*SIN (PHI) *SUMY-S*SUMZ) *ELMT 
CIRCR=(S*COS (PHI) *SUMX+S*SIN (PHI) *SUMY+C*SUMZ ) *ELMT 
CIRCP=(-SIN (PHI) *SUMX+COS (PHI) *SUMY ) *ELMT 

ENDIF 

RETURN 

END 


CR KKK XFUNCTION TAPER 


FUNCTION TAPER(RHO) 


C SPECIFY ANTENNA ILLUMINATION FUNCTION. REAL FUNCTION OF 
C RHO ONLY 


TAPER=1. 
RETURN 
END 
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